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Abstract 

Divide-and-conquer based methods for Bayesian inference provide a general approach for tractable 
posterior inference when the sample size is large. These methods first divide the data into smaller subsets 
and sample from the posterior distribution of parameters in parallel across all subsets, and then combine 
posterior samples from all the subsets to approximate the full data posterior distribution. Sampling in 
the first step is more efficient than sampling from the full data posterior due to smaller size of any sub¬ 
set. Since the combination step takes negligible time relative to sampling, posterior computations can 
be scaled to massive data by dividing the full data into sufficiently large number of data subsets. One 
such approach relies on the geometry of posterior distributions estimated across different subsets and 
combines them through their barycenter in a Wasserstein space of probability measures. We provide the¬ 
oretical guarantees on the accuracy of approximation that are applicable in many applications. We show 
that the geometric method approximates the full data posterior distribution better than its competitors 
across diverse simulations and reproduces known results when applied to a movie ratings database. 

Key words: barycenter; big data; distributed Bayesian computations; empirical measures; linear pro¬ 
gramming; optimal transportation; Wasserstein distance; Wasserstein space. 


1 Introduction 

Developing efficient sampling algorithms is an active area of research motivated by tractable Bayesian in¬ 
ference in large sample settings. Sampling remains a primary tool for inference in Bayesian models, with 
Markov chain Monte Carlo (MCMC) and sequential Monte Carlo (SMC) providing two broad classes of al¬ 
gorithms that are routinely used. Most MCMC and SMC algorithms face problems in scaling up to massive 
data settings due to memory and computational bottlenecks that arise; this has motivated a rich literature 
in recent years proposing a variety of strategies to enable better performance in such settings. Our focus is 
on proposing a very general divide-and-conquer technique, which is designed to combine results from any 
posterior sampling algorithm applied in parallel using subsets of the data. 

Massive data pose three major problems for existing sampling algorithms. First, if full data require 
multiple machines for storage, then any sampler can access only a small fraction of the data and posterior 
sampling given the full data is expensive due to extensive communication among machines. Second, if full 
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data are available to the sampler, then sampling is infeasible because computation of Hessians and accep¬ 
tance ratios often scale as 0(n^], where n is the sample size. Finally, sampling in hierarchical Bayesian 
models requires generation of 0(n) latent variables, which becomes inefficient as n increases. A variety of 
methods exist to address these issues using optimization and sampling. 

Optimization-based methods for Bayesian inference obtain an analytic approximation of the full data 


posterior distribution. The two most common techniques are polynomial approximation ( Rue et aLj 20091 
and projection of the full data posterior distribution on a class of distributions with analytically tractable 


posterior densities, which includes variational Bayes and expectation propagation (Wainwright and Jordan 


2008 Gelman et al.[|2014 l. Both techniques estimate parameters of the approximate distribution using a 
variety of optimization algorithms ( Thn and Nott} 2013 Wand 20141. Stochastic approximation significantly 
improves the efficiency of estimation by accessing the data in small batches and updating the parameter 


estimates sequentially (Broderick et al.[ 2013 Hoffman et al. 20131; however, optimization can be nontrivial 


for complex likelihoods frequently used in hierarchical models. Further, variational Bayes and expectation 
propagation often have excellent predictive performance but can be highly biased in estimation of posterior 
uncertainty and dependence. 

There is extensive work in sampling-based methods for Bayesian inference. The three main techniques 
used in these methods are as follows. First, subsampling-based methods obtain posterior samples condi¬ 
tioned on a small fraction of the data ( Maclaurin and Adams[ 2015| ). Coupling of subsampling with modified 
Hamiltonian or Langevin Dynamics improves posterior exploration and convergence to the stationary distri- 


bution (jWelling and Teh[ 201 1[ Ahn et al. 

20 12^ Korattikara et al.[|2014 Lan et al. 

2014 

Shahbaba et al.| 

2014 1 ; see Bardenet et al. 

(20151 for a review. Second, the exact transition kernel in posterior sampling is 


replaced by an approximation that significantly reduces the time required to finish an iterafion of the sampler 
(Johndrow et al. 2015| Alquier et al.[ 2016| ). Finally, divide-and-conquer approaches first divide the data 
into smaller subsets and sample in parallel across subsets, and then combine the posterior samples from all 
the subsets. Our focus is on scalable Bayesian methods based on the divide-and-conquer technique. These 
methods have two subgroups that differ mainly in their sampling scheme for every subset and their method 
for combining posterior samples obtained from all the subsets. 

The first subgroup modifies the prior to sample from the posterior distribution of the parameter condi¬ 
tioned on a data subset. Let k be the number of subsets, 7t(0) be the prior density of parameter 0, and li(0) 
be the likelihood for subset i (i = 1,..., k). Samples from subset posterior distribution i are obtained using 
li(0) and 71(0]^/^^ as the likelihood and prior. Consensus Monte Carlo combines subset posterior samples 
by averaging (Scott et al. 2016|). This relies heavily on the normality assumption, which is relaxed using a 


combination based on kernel density estimation ( [Neiswanger et al.[|2014[ ). Both methods perform poorly if 
the supports of subset posteriors are different, which motivates the combination using the Weierstrass trans¬ 
form and random partition trees ( |Wang and Dunson 2013[ Wang et al. 20151. These methods offer simple 
approaches to combine samples from subset posterior distributions but have two major limitations. First, the 
sampling algorithm depends on the model parameterization. Second, transforming 7r(0) to 7t(0)^/’^ makes 
posterior computations harder because existing off-the-shelf sampling algorithms cannot be used directly. 

The second subgroup modifies the subset likelihood to sample from a subset posterior distribution and 
combines samples from subset posterior distributions through their geometric center. These methods modify 
the likelihood to li(0)’^ and use prior 7t(0) to sample from subset posterior distribution r (i = 1,..., k). 
M-Posterior combines subset posterior distributions through their median in the Wasserstein space of order 


1 (Minsker et al. 20141. The robustness of the median implies that it could ignore valuable information 


in some subset posterior distributions, which motivates combination through the mean in the Wasserstein 
space of order 2 called Wasserstein Posterior (WASP) ( Srivastava et al.[[2015[ ). The WASP approach strikes 
a balance between the generality of sampling and the efficiency of optimization; however, its computations 
are developed for independent identically distributed (iid) data and its theoretical properties are unknown. 

Our main goal is to study three theoretical properties of WASP and apply WASP in a variety of practical 
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problems. The iid assumption of WASP rules out many important practical problems, including regression 
and classification, where the data are independent and non-identically distributed (inid). We relax this 
assumption and our results are applicable to any inid data. Second, we show that if the number of subsets 
are chosen appropriately, then the WASP achieves almost the same rate of convergence as that of the full 
data posterior distribution. This implies that WASP can be used as an efficient alternative to the full data 
posterior distribution for uncertainty quantification in massive data settings. Third, we show that the method 
for estimating WASP is independent of the form of the model, so WASP is very general and can be easily 
used for estimating posterior summaries for any functional of the model parameters. We emphasize that 
WASP is not a new sampling algorithm but a general approach to easily extend existing sampling algorithms 
for massive data applications. 


2 Preliminaries 

2.1 Wasserstein space, Wasserstein distance, and Wasserstein bary center 

We recall elementary properties and definitions related to the Wasserstein space of probability measures. 
Let (©, p) be a complete separable metric space and T(0) be the space of all probability measures on ©. 
The Wasserstein space of order 2 is defined as 


T 2 (©] := 


|qGT(©]: 


p(0o, 0]^p.(d0) < oo 

0 


( 1 ) 


where 0o G © is arbifrary and T’ 2 (©) does nof depend on fhe choice of 0o. The space T 2 (©) is equipped 
wifh a nafural disfance befween ifs elemenfs. Lef p, "v G T 2 (©) and n(p, y) be fhe sef of all probabilify 
measures on © x © wifh marginals p and y, fhen fhe Wassersfein disfance of order 2 befween p and y is 
defined as 


W2(p,y) 


inf 

7ter[(|j,,v) 


0X0 


d^(x,y) d7r(x,y] 


( 2 ) 


In our applicafions p is fhe Euclidean mefric and we refer fo T 2 (©) and W 2 as fhe Wassersfein space and fhe 
Wassersfein disfance wifhouf explicifly mentioning fheir order. If fTi,..., Flic are a collecfion of probabilify 
measures in T 2 (©)^ then fheir barycenfer in T 2 (©) is defined as 


n = 


argmin Y_ 


(3) 


ney2(0] 


j=i 


This generalizes fhe Euclidean barycenfer, which is fhe sample mean, fo T 2 ( 0 ] (Agueh and Carliei] 2011 1 . 


The barycenfer fT is analyfically infracfable, excepf in few special cases. Eel 6 q(x) = 1 if a = x and 


0 ofherwise. If Xi 


ji,... , Xjm, are samples from fTj (j = 1 ,..., k), fhen flj (•] = 6 xjt(-]/Tn- is an 

empirical measure fhaf approximafes fTj (j = 1,..., k). If fl is assumed fo be an empirical measure, fhen 
fhe opfimizalion problem in Q reduces fo a linear program; see Cufuri and Doucef] ( |2014 ), Cartier el al. 
(20151, and Srivasfava ef al.|([2015[) for differenl algorilhms fo solve Ihis linear program. 


2.2 Stochastic approximation and subset posterior density 

Consider a general sel-up for inid dala. Eel = (Yi, ■ • •5 Yn) be n observalions and fhe dislribulion 
of Yt is Pe t, i = 1,..., n, where 0 ties in fhe parameter space © C Assume lhal Pe^i has density 
pt(-!0) wifh respecl fo fhe Eebesgue measure, so dP 0 ^t(yt) = pt(yi|0)dyt and fhe likelihood given Y^^^ 
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is 1(0) =Or=iPi(yil®^ Given a prior distribution TT on 0 that has density tt with respect to the Lebesgue 
measure, the posterior density of 0 given using Bayes theorem is 


^(9iyM)^ , nr.iPiitiii8W) 


i(e)ji(0) 


Jo nr=i Pi(yi I 0)7r(0)d0 Jq l(0)7t(0)d0' 


(4) 


In most cases 7t(0 | is analytically intractable, and accurate approximations of 7 t( 0 | Y^^^) are 

obtained using Monte Carlo methods, such as importance sampling and MCMC, and deterministic approxi¬ 
mations, such as Laplace’s method and variational Bayes. For example, in the context of logistic regression, 
Pe i is the Bernoulli distribution with mean 1/ |l -|- exp(—x40)|, where xj is the ith row of the design 
matrix X G and 0 = The posterior density of 0 is analytically intractable, and it is typical to rely 

on Gibbs samplers based on data augmentation. These samplers introduce latent variables {zi, 1 = 1,..., n} 
and alternately sample the latent variables and the parameters from their full conditional distributions. Re¬ 
lated algorithms are very common and are computationally prohibitive for large n because they require 
repeated passes through the whole data. 

Divide-and-conquer-type methods resolve this problem by partitioning the data into smaller subsets. 
Let k be the number of subsets. The default strategy is to randomly allocate samples to subsets. Let 
Y[j] = = (Yji,..., Yjraj) denote data on the jth subset, where mj is the size of the jth subset 


and = TT- We assume that mj = m (j = 1,..., k) for ease of presentation, so n = km, the 

likelihood given Y[j] is lj(0) = OrJi PjilyjilQ)^ and 1(0) in Q equals nf=i 
and Srivastava et al. (20151 define subset posterior density j given Y[j] as 


Minsker et al. 


(20141 


^ .QIY . {n[^iPji(yjil9)iM9) ij(9)Me) 

^ j0{niiiPji(yji|0)F7r(0)d0 j0lj(0)Y7r(0)d0’ 


where y is a positive real number such that giym ^ n ^ 927^1 for some gi, 92 > 0- In the present 
context, we assume that y = k with gi = 92 = 1 following [Minsker et al. (20141; more general conditions 
on y are defined lafer in Secfion 3.2 This modified form of subsef posferior compensafes for fhe facf fhaf jfh 
subsef has access fo only (m/n)-fracfion of fhe full dafa and ensures fhaf 7trn,(0 I Y[j]) and 7 rTt (0 | Y^^^) in 
(Q have variances of fhe same order. Minsker ef al.| ( |201^ refer fo fhis as stochastic approximation because 
raising Ij ( 0 ) (j = 1 ,..., k) fo fhe power y is equivalenf fo replicafing every Xji (i = 1 ,..., m) y-fimes so 
fhaf 7rTa(0 I Y[j]) (j = 1,..., k) are noisy approximafions of 7 r (0 | Y^^^). 

One advanfage of using stochastic approximation fo define 71^(0 I Y[j]) in ([^ is fhaf off-fhe-shelf 
sampling algorifhms can be used direcfly even when fhe prior density is fhe form of a discrefe mixture. 
Consider a simple example of univariate density estimation using Dirichlet process mixtures of Gaussians. 
Let Xt (i = 1,..., n) be iid samples from a distribution Pq with density po. The data are randomly split into 
k subsets of equal size m. The truncated stick-breaking representation of Dirichlet process prior implies that 
the prior distribution FI on T has a finite mixture representation, where T is the set of probability distributions 
that have a density. The subset posterior density j of our competitors is proportional to Ij (p)7r(p)^/'^, where 
p is the density of a distribution P G IP, lj(p) is the likelihood for subset j, and 7r(p) has the form of 
a Dirichlet process mixture. In this case, existing Gibbs sampling algorithms cannot be applied directly. 
In contrast, we show in Appendix that modification of the likelihood using stochastic approximation in 
WASP allows the use of existing Gibbs sampling algorithm for density estimation. 

Stochastic approximation does not add any extra burden to the computations required for sampling 
from the subset posterior distribution of 0 conditioned on m observations. We raise the likelihood in every 
subset to the power y. This is equivalent to replicating observations y-times, which seems to offset the 
benefits of partitioning. However, the replication of observation is not required in implementation of the 
sampler; we simply modify the likelihood in the full data sampler by raising it to the power y. For example. 
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Figure 1: Binned kernel density estimates of full data posterior distribution, subset posterior distributions, 
and WASP for coefficients (0i, 62 ) in logistic regression. The x and y axes represent posterior samples for 
01 and 02 - The true values of 0i and 02 are —1 and 1 (black triangle). 


stochastic approximation is easily implemented using the increment_log_prob function in Stan (Stan 
Development Team| 2014| ). We provide more examples for a variety of models in Section]^ 

A simple logistic regression example demonstrates that 7rTa(0 I Yyj) in ([^ is a noisy approximation of 
7 t (0 I ) in (0. We simulated data for logistic regression with n = 10^, p = 2, 0 = (—1,1)^, and entries 
of X randomly set to ±1 (Figure[^. We set y = k = 40 and obtained samples of 0 from 7 t (0 | and 
from 7tTn(0 I Y[j]) (j = 1 ,..., k) using the Stan’s HMC sampling algorithm. The contours for the subset 
and full data posterior densities are very similar, indicating all densities have similar spreads. We also notice 
that subset posteriors are noisy approximations of the full data posterior in that most of them have a bias and 
do not concentrate at the true 0 . 


3 Wasserstein Posterior (WASP): The general framework 


3.1 Definition and estimation of the WASP 

The WASP approach combines subset posterior distributions fTrat- I Yyj) (j = 1 ,... ,k) through their 
barycenter in T’ 2 ( 0 ]> where the density of nTa(- I Yyj) is 7rTa(- I Y[j]) in ([^. The barycenter represents a 
geometric center of a collection of probability distributions that can be efficiently computed using a linear 
program. Motivated by this, Srivastava et al.j ( 2015| ) proposed to combine a collection of subset posterior 
distributions through their barycenter in the Wasserstein space called WASP. Assuming that subset posterior 
distributions nTn(- I Y[j]] (j = 1 ,..., k) have finite second moments, the WASP is defined using (|^ as 


1 , 


nn(- I Y(-5] = argmin Y_ rW2{n,n^(- | Y[p)}. 


( 6 ) 


neyato] 


i=i 


The WASP is analytically tractable only in special cases, but it can be estimated using a linear program 


if the subset posterior distributions have an atomic form (Srivastava et al. 20151. Let {0ji,..., 0js} be 
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the 0 samples obtained from subset posterior density j in Q using a sampling algorithm, including HMC, 
MCMC, SMC, or importance sampling. Approximate jth subset posterior distribution TTml' I Yyj) using 
the empirical measure 


s 1 


i=l 


Srivastava et al. (20151 approximate the WASP as 

k S 

j=i i=i 


0 ^ aji ^ 1, 


k S 

^ ^ O-ji = 1) 
j=li=l 


(V) 


( 8 ) 


where aji (j = 1,..., k; i = 1,..., S) are unknown weights of the atoms. There are many specialized 
algorithms to estimate the WASP that exploit the structure of the linear program in ([^ when nTn,(- I Y[p) 
and TTn(- I are restricted to have atomic forms in (|7]l and respectively; for example. 


Cuturi 


and Doucet (2014) extend the Sinkhorn algorithm of Cuturi|([2013 1 using entropy-smoothed sub-gradient 


methods, Carlier et al. (2015) develop a non-smooth optimization algorithm, and Srivastava et al. (2015) 
propose an efficient linear program that exploits the sparsity of constraints to solve (|^. 


3.2 Theoretical properties of the WASP 


We study asymptotic properties of the WASP defined in ([^ thaf are important for applications in scalable 
Bayesian inference. First, the rate of concentration of the WASP around the true value of the parameter is 
unknown. Let Oq be the true parameter, then the full data posterior distribution Fin (• | ) converges to 60 ^ 

in Hellinger distance at the optimal parametric rate of up to logarithm factors (Ghosal et al. 2000); 

however, similar guarantees are unknown for the WASP. In addition, in many applications, the parameter 0 
is not of immediate interest and the interest lies in f (0) for some known function f. We do not know if the 
WASP of k subset posterior distributions for f(0] converges to 6 f ( 0 ^). The current asymptotic justification of 


WASP relies only on posterior consistency under an iid assumption; for example. Theorem 3.3 in Srivastava 


et al. (2015|l. We obtain much stronger asymptotic support by showing rates in an inid case, including for 


functionals of the original parameters. 

Our theoretical setup is based on the following assumptions. 

(Al) 0 is a compact space in p metric, and 0 o is an interior point of 0 ; giym ^ n ^ g 2 yTU for some 
constants gi, 92 > 0 . 


(A2) For any 0,0 G 0 and j = 1, • 


., m, there exist positive constants a and Cl, such that 

T0,0')^CLp2“(0,0'), 




where (0,0 ) is the pseudo Hellinger distance defined in 


A.l 


(A3) (Entropy Condition) There exist constants Di > 0 , 0 < D 2 < Df/ 2 ^ 2 , a function ¥(u,r) ^ 0 that 
is nonincreasing in u G and nondecreasing in r G M^, such that for all j = 1,... ,k, for any 
u, r > 0 and for all sufficiently large m. 


H[] (u, {pj(y|0) : 0 G 0,h.Tai(0,0o) ^ t} ,h.Tnj) ^ 'T(u,r] for all j = 1,... ,k; 


and 


Dir 


Dit2/212 


a/¥(u, r)du < D 2 


mr 


where pj (y |0) = {pji(yji I 0], • ■ • iPjmjyjra I 9)}^ and H[] is the hm.) -bracketing entropy of the set 
{pj(y!0) : 0 G 0 , Hraj (0, 0o) ^ f}, which is defined in[ 


A.2 
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(A4) (Prior Thickness) There exist positive constants k and c^, such that uniformly over all j = 1,..., k, 


n 




1 = 1 


, Pj^(Yji|0o)^ , ^ log' 



exp(—Cyrklog^ m] 


where log_|_ x = max(logx, 0) for x > 0. 


(A5) The metric p satisfies that for any N G N, 0i,..., 0 n , 0^ £ 0 and nonnegative weights Wi = 
1 , 


N 


N 


P ^Wi0i,0' ^ y Wip(0i,0O 


a=i 


i=l 


For theory development, we have assumed compact support in (Al) and lower bounded pseudo Hellinger 
distance in (A2), which is similar to Theorem 10 in Ghosal and van Per Vaai^ (2007 1. Typically, a = 1 
for most regular models; e.g., generalized linear models. If the model is non-regular, then a can be less 
than 1. For example, the densities may have discontinuities depending on the parameter (see [Ibragimov and 


Has’minskii (19811 Chapter V and VI). (A3) parallels the entropy condition used in Theorem 1 of Wong 
and Shen] ( 19951, which has been adapted here for the inid setup using the generalized bracketing entropy, 
and will simplify to a similar entropy condition to that in Theorem 1 of Wong and Shen] ( 19951 if the data 
are iid. The convexity property of p in (A5) is mainly used to establish an averaging inequality under W 2 
distance and is satisfied by for example, fhe Euclidean mefric and Iq mefric wifh q ^ 1; see Lemma [ET] in 
fhe Supplemenfary Material. 

The fwo fheorems below describe fwo asympfofic properties of WASP for models fhaf safisfy assump¬ 
tions (A1)-(A4) above. The asympfofic behavior of fhe WASP for inid dafa is described in fwo sfeps. The 
firsf fheorem describes fhe asympfofic behavior of subsef posteriors disfribufions nm,)- | Y[j]) (j = 1,..., k). 
The second fheorem describes fhe asympfofic behavior of fhe WASP using fhe asympfofic behavior of subsef 
posterior disfribufions. 


Theorem 3.1 If Assumptions (A1)-(A4) hold for the ]th subset posterior nnx(' I Y[j]) with ] 
then there exists constants a and Ci that do not depend on j, such that as m ^ 00 , 


= 1 , 


,19 


E[Wi{n^(- I Y[q],6eo(-)}] ^ Cl 


log^ m 


m 


(9) 


where expectation is with respect to probability measure Pg^ 


m) 


Theorem 


3.1 


proves posferior convergence in expecfafion Epg^, which is sfronger fhan fhe commonly 
sfudied posterior convergence in probabilify. Assumpfion (A4) is crucial in providing a sfronger confrol 
over fhe fail probabilify as fhe posferior probabilify mass moves away from fhe frue parameter 0o, fypically 


wifh an exponentially decaying rale. For regular models wifh a = 1, Theorem 3.1 says fhaf fhe W 2 dislance 
belween fhe subsef posferior Flral- I Y[j]) and fhe delfa measure af 0o shrinks lo zero al a nearly paramefric 
rate up lo some logarilhm faclors. The rate is fhe opfimal paramefric rale because each 

subsef only has sample size m. For nonregular models wifh a < 1, fhe convergence rale can be faster 
fhan paramefric, which achieves fhe same rale as fhe frequenfisl maximum likelihood eslimalor up lo some 


logarilhm faclors; see Ghosal el al. (19951, Ibragimov and Has’minskii (19811. A combination of k disjoinl 
subsels yields fhe following convergence resull for fhe WASP. 
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Theorem 3.2 If Assumptions (A1)-(A5) hold for all subset posteriors \ Y[j]) with] = then 

<35 m —7- oo, 


W2{nn(-1 



( 10 ) 


where Op is in -probability. 


Theorem 3.2 shows that the WASP of k subset posterior distributions converges to the true parameter Oq at 
a rate of m“2« up to logarithm factors. When a = 1, as in regular models, this rate is about the same as 
the rate obtained in Corollary 3.12 and Theorem 3.14 of [Minsker et al. ( 2014| l for the M-Posterior. If a = 1 
and the number of subsets k increases slowly with n, for example k = O(log‘^ n) for some constant c > 0, 


then Theorem 


3.2 


implies that the WASP has a near optimal convergence rate Or 


optimal parametric rate on the full dataset is Op (n . 


log‘= 


since the 


Suppose f : 0 I— )■ M*! is a function that maps 0 to {fi(0], ■ ■ •, fq (0)}, where q ^ 1 is a positive integer. 


A direct application of Lemma 8.5 in Bickel and Freedman (19811 gives the following corollary about the 
WASP of a function of 0. As long as the function is bounded almost linearly by the p metric in ([T]), its 
WASP possesses the same posterior convergence rate as in Theorem|3.2| 


Corollary 3.3 Suppose f(-) = {fi(0; ■ ■ • ; fq (■)} a function that maps 0 i—)> such that |f(0)P = 
^?^^{ft(0)}^ ^ Cf{l + p^(0, 0o)}. where Cf > 0 is a fixed constant. If the conditions in Theorem 
and fttflTT^)- I Y^^^) represents the WASP of the subset posterior distributions for f (0), then a5 m —)• oo, 


hold 


W2 [m. 


Y( 


TL h C 


(.)} 


= 0 , 


'log 


2/a 


m 


m 


l/a 


3.3 Computation of the WASP 


Corollary 3.3 is very useful in applications because it says that the combination step in the WASP is 
independent of the model parametrization. Let fttnrn,(- I Y[j]) be the jth subset posterior distribution for f (0) 
(j = 1,..., k), then the WASP of k subset posterior distributions converges to f (0o) at the rate obtained in 


Theorem 3.2 In practice, we have Sj posterior samples of 0 obtained from subset posterior j denoted as 0j| 
(i = 1,..., Sj; j = 1,..., k). Algorithm [^estimates an atomic approximation of ft)fTTT^(- | YI"-!), denoted 

asfUn^(- I YI^I), based on the subset posterior samples f(0ji) (i = = l,...,k). The atomic 

form of the WASP is supported on a grid with mesh-size e estimated from the subset posterior samples of 
f (0). Algorithm [^estimates the weights of the atoms located on the grid by solving a discrete version of Q. 
The support of ftjnT^(- | Y1"-1) scales exponentially in k, but theoretical properties of discr ete barycenters 


imply that f (jn,- 


'n V 
Y(i^) 


is supported only on 0(k) elements of the grid; see Theorem 2 in 


Anderes et al. 


(20161. We exploit this sparsity by adapting the algorithm in Srivastava et al. (20151 and by relying on 


Gurobi to exploit sparsity in linear programs ( |Gurobi Optimization Inc. 20141. A general algorithm that 
exploits the sparsity in the support of the atomic approximation of WASP is left for future research. In all 
our experiments, ftin^(- | Y^^^) provides an approximation of the full data posterior distribution of f(0). 



































Algorithm 1 Estimation of the WASP for f (0) in Corollary |3.3|given samples of 0 from k subset posteriors 


Input: Samples from k subset posteriors, {0ji : Gjr ~ | Yy]), i = 1,..., Sj, j = 1,. .., k}; mesh size e > 0. 

Do: 

1. Define 4)1 = (4){]^, ■ • •, 4>{q) = f(0jll P “ j = the matrix of atoms of subset posterior j, <Dj G 

with cf)? as row i (i ^ 1, . . . , Sj). For T ^ 1, . . . , q. let (j>inin ^ (cl>min l, ■ • ■ , 4>min q ] with cj^min r ^ min and 4)niax ^ 


( 4^ max 1) • • ■ ! 4^ max q } with c|) max r — Hiax ^ . 


2. Set the number of atoms in the empirical approximation for the WASP g = giX...xgq, where g r = 


_ r tt’max r 4>miii r 1 /•-*- _ 




3. Define the matrix of WASP atoms O G ^ with rows formed by stacking vectors 

4^ min 1 H“ ( 4^niax 1 4^min 1 } i • • • ) 4^ min q H“ (4^ max q 4^ min q)^5 (^r — 


4. Set the distance matrix between the atoms of WASP and the jth subset posterior, Dj G , as 

q _ 

(Dj)uv = (u. = I,- - •, g; V = 1,... ,Sj; i = 1,... ,k), 

r=l 

where 4 >ut is the (u, r)-entry of O. 

5. Estimate dp,..., dg by solving the linear program iH} in Appendix[D| 

Return: ftin(- | atomic approximation of fJEl^ (■ | ]. 


4 Experiments 

4.1 Setup 


We compared WASP with consensus Monte Carlo (CMC) (Scott et ^ 20161, semiparametric density prod¬ 


uct (SDP) ( [Neiswanger et ah' 20141, and variational Bayes (VB). The sample sizes and the number of 
parameters in our experiments were chosen such that sampling from the full data posterior distribution was 
computationally feasible. Every sampling algorithm was run 10,000 iterations. We discarded the first 5,000 
samples as bum-in and thinned the chain by collecting every fifth sample. Convergence of the chains to 
their stationary distributions was confirmed using trace plots. All experiments were run on an Oracle Grid 
Engine cluster with 2.6GHz 16 core compute nodes. Eull data posterior computations were allotted memory 
resources of 64GB, and all other methods were allotted memory resources of 16GB. 

The sampling algorithm for the full data posterior was modified fo obtain samples from the subset 
posteriors in CMC, SDP, and WASP. The sampling algorithms for subset posteriors in CMC and SDP were 
the same and were based on Equation (2) in Scott et al. ( |2016 ). The sampling algorithm for subset posteriors 
in WASP was based on Q. Samples from the approximate posterior distributions of 0 in CMC, SDP, and 
WASP were obtained in two steps. Eirst, samples from subset posteriors of 0 were obtained in parallel across 
k subsets. Second, the samples of 0 from all the subsets were combined using implementations of CMC and 
SDP in parallelMCMC package (Miroshnikov and Conlon 2014| and using Algorithm [T] for the WASP. 

The full data posterior distribution obtained using MCMC served as the benchmark in all our compar¬ 
isons. Eet 7t(0 I be the density of the full data posterior distribution for 0 estimated using sampling 

and 7t(0 I be the density of an approximate posterior distribution for 0 estimated using the WASP or 
its competitors. We used the following metric based on the total variation distance to compare the accuracy 
7t(0 I Yt^^) in approximating 7 t(0 | Y^^^) 


accuracy 1^(0 


/(n) 


>} 


= 1 - - 

2 


0 


7t(0|Y'^5)-7r(0|Yt’^’; 


d0. 


( 11 ) 
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The accuracy metric lies in [0,1] (Faes et al. 20121. The approximation of full data posterior density by ft is 


poor or excellent if the accuracy metric is close to 0 or 1, respectively. In our experiments, we computed the 


kernel density estimates of ft and n from the posterior samples of 0 using R package KernSmooth (Wand 
2015 1 and calculated the integral in ( [TT] ) using numerical approximation. 


4.2 Simulated data: finite mixture of Gaussians 

Finite mixture of Gaussians are widely used for model-based classification, clustering, and density estima¬ 
tion ( Fraley and Rafter^ 20021. Let n, p, and L be the sample size, the dimension of observations, and the 
number of mixture components. If iJi G is the rth observation (r = 1,... ,n.), then the mixture of L 
Gaussians assumes that any y G ..., is generated from the density 


fmix (y I 0) = ^ 7T; Np (y I (X;, 10, 


( 12 ) 


t=i 


where n = (rti,..., ttl) lies in a (L — l]-simplex, p,^ and (I = 1,..., L) are the mean and covariance 
parameters of a p-variate Gaussian distribution, and 0 = {n, ..., (Xp, Zi,..., Zl}. We set L = 2 and 

p = 2 and simulated data from (|43|) using 7t = (0.3,0.7], (Xi = (1,2)^, p ,2 = (7,8]^, and Z; = Z 


(1 = 1, 2), where Z 12 = 0.5, Zn = 1, and Z 22 = 2. We performed 10 simulation replications. 

This simple example demonstrated the generality of WASP in estimating the posterior distribution of 
functions of 0 as described in Corollary|3.3[ We defined fwo nonlinear funcfions of 0 as 


Pt = (It)l2/{(It)ll(^t)22}^^^ 1 = 1,2, g(x) =fi-nix{(x,x)'} XG 


(13) 


where pt is the correlation of Ith mixture component and g(x] is the value of density fmix in ( |43l ) when 
y = (x,x)^. Our simulation setup implied that pi = p 2 and g(x) was bimodal for x G M. We completed 
the hierarchical model in ( |43| ) by specifying independent conjugate priors on n and (p.;, Zj,) (I = 1, 2) as 

7t ~ Dirichlet(l/2,1/2), p.^^ | Zi^ ~ 3^2(0, lOOZ^), Z; ~ Inverse-Wishart(2, 4 I 2 ) (1 = 1,2), (14) 

where 2 is the prior degrees of freedom and 41p is the scale matrix of the Inverse-Wishart distribution. 


The posterior samples of 0 were obtained using Gibbs sampling (Bishop 20061, which were used to obtain 
posterior samples for pi, p 2 , and g. 

We compared WASP with the posterior distributions estimated using CMC, Gibbs sampling, SDP, and 


VB. We used the VB algorithm developed in Bishop (20061. Two values of k G {5,10} were used for 


CMC, SDP, and WASP and full data were partitioned into k subsets such that the mixture proportions 
were preserved in every subset. The approximate posterior distributions of pi, p 2 , and g(x),x G M, under 
each method were estimated using the subset posterior samples obtained after modifying the original Gibbs 
sampler. The sampling algorithm for WASP is described in Section [2T] of Supplementary Material. 

We compared the accuracy ( [TT] ) of CMC, SDP, VB, and WASP in approximating the full data posterior 
distributions of pi, P 2 , and point-wise 90% credible bands of g(x) for x G M, denoted as go.o 5 (^) and 
90.95(x)- CMC, SDP, and WASP accurately approximated the full data posterior distributions of pi and P 2 
for both ks, but VB underestimated the posterior uncertainty in pi and p 2 . CMC, VB, and WASP were 
very accurate in estimating go.o 5 (^) and go.gs)^) for x G M, whereas the application of SDP failed due to 
a numerical error in matrix inversion (Table [TJ. This provides an empirical verification of Corollary |3.3[ 
showing that the accuracy of the WASP was unaffected by the form of the parameters in the combination 
step. Theoretical guarantees similar to Corollary |3.3| were unavailable for CMC or SDP, but our numerical 
results illustrated that a similar result might also hold for these methods in mixture models. 
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Table 1: Accuracies of the approximate posteriors for pi, p 2 , and go.o 5 (^] and go. 95 (x) for x G M. The 
accuracies are averaged over 10 simulation replications. Monte Carlo errors are in parenthesis. CMC, 
consensus Monte Carlo; SDP, semiparametric density product; VB, variational Bayes; WASP, Wasserstein 
posterior 



Pi 


P2 

90.05 

90.95 

VB 

0.77 (0.31) 

0.76 (0.29) 

0.99 (0.00) 

0.99 (0.00) 


k = 5 

k = 10 

k = 5 

k= 10 

k = 5 

k= 10 

k = 5 

k = 10 

CMC 

0.97 (0.01) 

0.96 (0.01) 

0.96 (0.01) 

0.96 (0.01) 

0.99 (0.00) 

0.99 (0.00) 

0.99 (0.00) 

0.99 (0.00) 

SDP 

0.97 (0.01) 

0.96 (0.01) 

0.95 (0.01) 

0.96 (0.01) 

- 

- 

- 

WASP 

0.97 (0.01) 

0.95 (0.01) 

0.97 (0.01) 

0.96 (0.01) 

0.99 (0.00) 

0.99 (0.00) 

0.99 (0.00) 

0.99 (0.00) 


4.3 Simulated data: Linear mixed effects model 


Linear mixed effects models are extensively used in extending linear regression to account for longitudinal 
and nested dependence structures. Let n, s, and st be the sample size, total number of observations, and 
total number of observations for sample i (i = 1,..., n) so that s = Si. Suppose Xt G RSixP and 

Zi G include predictors in the fixed and random effects components, respectively. Letting TJi G 

be the response for sample i, the linear mixed effects model assumes that 


iJt I P,Ui,T2 ~ AfsilXi |3+ZiUi,T2lT 


Ui 




fi = L 


,rL 


(15) 


where Ui G is the random effect for sample i with mean 0 and r x r covariance 1, p G denotes the 
fixed effecfs, and is fhe error variance. The model paramefers are 0 = {p, L, t^}. 

We simulated dafa for a fixed n and s and varying p and r. We chose fwo values of (p,r) G 
{(4,3), (80,6)}, fixed n and s fo be 5000 and 100,000, and randomly assigned fhe s observafions fo n 
samples. The fwo choices of (p,r) ensured fhaf fhe number of unknown paramefers in P and L was 10 
and 100 in fhe former and laffer cases. The enfries of Xt and Zi were sef fo 1 or —1 with equal probabil¬ 
ity for every i. We fixed P enfries as —2 and 2 alfernafely, = 1, and L = LL^, where L was a lower 
friangular mafrix wifh Lti = -v/i (1 = 1,..., r) and fhe off-diagonal enfries from L 21 fo LT-(T._i) increasing 
linearly from \/r~^/100 fo Y^(r — 1)/100. We used fhis sefup fo simulate dafa from ( fT5| ) and performed 10 
replicafions. 

We used the HMC algorithm in Stan for sampling from the full data and subset posterior distributions. 
The full data posterior computations were feasible for the chosen values of n and s and posterior samples 
were obtained after completing the hierarchical model in ( fT5| ) by using the default weakly informative priors 
for p, L, and in Stan. Two values of k G {10, 20} were used for CMC, SDP, and WASP, and the n samples 
were randomly partitioned into k subsets. The sampling algorithms for subset posterior distributions for the 
three methods were implemented in Stan and posterior samples of 0 were obtained in parallel across k 
subsets. This was followed by a combination step to estimate the approximate posterior distributions for the 
three methods. The sampling algorithm for WASP is described in Sectionof Supplementary Material. 

We compared the accuracy ( [TT] ) of CMC, SDP, VB, and WASP in approximating the marginal posterior 
distributions of fixed effecfs, variances and covariances of random effecfs, and the joint posterior distribu¬ 
tions of three pairs of covariances of random effects. We used the streamlined algorithm for estimating the 
VB posterior for (3 and L ( Lee and Wand[ 20161. The four methods were significantly faster than the full 
data posterior distribution, while using only 25% of the memory resources, with VB being the fastest. CMC, 
SDP, VB, and WASP provided accurate approximations of the marginal posterior distributions of fixed ef¬ 
fecfs and covariances of random effecfs. The accuracy of CMC and SDP in approximating the marginal 
posterior distributions of variances of random effects depended on k and r, and VB underestimated the pos¬ 
terior variances of random effects. In these cases, the accuracy of WASP was stable for every k and r (Tables 
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[^and[^. All four methods showed similar accuracies in approximating the true joint posterior distributions 
of three pairs of covariances of random effects. The differences in accuracies for different values of k and 
r were due to the differences in numerical approximation of ( [TT| ) using kernel density estimation (Tables 
andj^and Figures]^ and |^; see Table [T0| in the Supplementary Material. 

The accuracy of CMC, SDP, and WASP decreased when k increased from 10 to 20 because subset pos¬ 
terior distributions conditioned on a smaller fraction of the data. This provided an empirical verification of 
Theorems 3.1 and 3.2 for the WASP. Our numerical results illustrated that a similar result might also hold 
for CMC and SDP. The stable performance of WASP compared to that of CMC and SDP in the approxi¬ 
mation of the posterior distributions of variances of random effects showed that the validity of the normal 
approximation for subset posterior distributions was crucial in obtaining accurate approximations of full 
data posterior using CMC and SDP. On the other hand, WASP results were free of any such assumptions 
and were valid for any nonlinear function of p, and Z; see Corollary|3.3| 


Table 2: Accuracies of the approximate posteriors for variances in ( [T5] ). The accuracies are averaged over 
10 simulation replications and across all diagonal elements of 1. Monte Carlo errors are in parenthesis. VB, 
variational Bayes; CMC, consensus Monte Carlo; SDP, semiparametric density product; WASP, Wasserstein 
posterior 



r = 

3 

r : 

= 6 

VB 

0.39 (0.21) 

0.48 (0.24) 


k= 10 

o 

csi 

II 

k = 10 

k = 20 

CMC 

0.93 (0.04) 

0.90 (0.06) 

0.85 (0.06) 

0.74 (0.1) 

SDP 

0.92 (0.05) 

0.83 (0.11) 

0.86 (0.08) 

0.77 (0.15) 

VB 

0.39 (0.21) 

0.39 (0.21) 

0.48 (0.24) 

0.48 (0.24) 

WASP 

0.97 (0.01) 

0.97 (0.02) 

0.97 (0.01) 

0.96 (0.01) 


Table 3: Accuracies of the approximate posteriors for covariances in ( [T5] ). The accuracies are averaged over 
10 simulation replications and across all off-diagonal elements of L. Monte Carlo errors are in parenthe¬ 
sis. VB, variational Bayes; CMC, consensus Monte Carlo; SDP, semiparametric density product; WASP, 
Wasserstein posterior 



r 

= 3 

r 

= 6 

VB 

0.95 (0.01) 

0.95 (0.02) 


k= 10 

k = 20 

k= 10 

k = 20 

CMC 

SDP 

WASP 

0.95 (0.01) 
0.93 (0.04) 
0.97 (0.01) 

0.95 (0.02) 
0.91 (0.06) 
0.96 (0.02) 

0.95 (0.02) 
0.91 (0.05) 
0.97 (0.01) 

0.94 (0.03) 
0.88 (0.08) 
0.96 (0.01) 


Table 4: Accuracies of the approximate two-dimensional joint posteriors for the covariances of random 
effects when r = 3 in ( [T5| ). The accuracies are averaged over 10 simulation replications. Monte Carlo errors 
are in parenthesis. CMC, consensus Monte Carlo; SDP, semiparametric density product; VB, variational 
Bayes; WASP, Wasserstein posterior 



(o'i2, Cia) 

(cTi2, 0'23) 

(cTia, 0'32) 

VB 

0.94 (0.01) 

0.94 (0.02) 

0.94 (0.02) 


k=10 k = 20 

k = 10 k = 20 

k = 10 

O 

CM 

II 

CMC 

SDP 

WASP 

0.93 (0.03) 0.93 (0.04) 

0.93 (0.03) 0.93 (0.04) 

0.95 (0.01) 0.95 (0.01) 

0.94 (0.02) 0.92 (0.03) 
0.92 (0.03) 0.91 (0.06) 
0.95(0.01) 0.94(0.01) 

0.93 (0.02) 
0.92 (0.03) 
0.95 (0.01) 

0.92 (0.03) 
0.90 (0.04) 
0.95 (0.01) 
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Table 5: Accuracies of the approximate two-dimensional joint posteriors for the covariances of random 
effects when r = 6 in ( fT5| ). The accuracies are averaged over 10 simulation replications. Monte Carlo errors 
are in parenthesis. CMC, consensus Monte Carlo; SDP, semiparametric density product; VB, variational 
Bayes; WASP, Wasserstein posterior 



(0'l2, CTia) 

(cTi2, 0'23) 

(O'la, 0'32) 

VB 

0.92 (0.01) 

0.94 (0.01) 

0.93 (0.01) 


k=10 k = 20 

k = 10 k = 20 

k = 10 

O 

(M 

II 

CMC 

SDP 

WASP 

0.93 (0.03) 0.91 (0.03) 

0.93 (0.03) 0.90 (0.05) 

0.94 (0.02) 0.94 (0.01) 

0.93 (0.02) 0.92 (0.02) 
0.93 (0.02) 0.90 (0.06) 
0.94(0.01) 0.94(0.02) 

0.93 (0.02) 
0.92 (0.03) 
0.94 (0.01) 

0.94 (0.02) 
0.92 (0.04) 
0.94 (0.02) 
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Figure 2: Kernel density estimates of the posterior densities of three covariance pairs when r = 3 in ( [T5] ), 
where Cab, acd on every panel represents the two-dimensional posterior density of (o^ab, o'cd)- MCMC, 
Markov chain Monte Carlo; CMC, consensus Monte Carlo; SDP, semiparametric density product; VB, 
variational Bayes; WASP, Wasserstein posterior. 


4.4 Simulated data analysis: Probablistic parafac model 

We use probabilistic parafac model as a representative example for nonparametric density estimation using 
WASP. Probabilistic parafac is an approach for nonparametric Bayes modeling of joint dependence in multi¬ 
variate categorical data (Dunson and Xing 20091. Let Xt = (xn,..., Xij,..., Xip ] be the data from sample 
i, where Xij has dj possible categorical values in {1,..., dj} (j = 1,... ,p). The hierarchical model for xij 
(i = j = l,...,p)is 


Xij 




H=1 


OO OO 

Zi - ^ Vh F[(l = L ~Beta(l,a), 


(j) 


.(j) 


h=l 


1<H 


H=1 


- Dirichlet(aji,..., ajdj], a ~ Gamma(aoc, ba), 


(16) 
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Figure 3: Kernel density estimates of the posterior densities of three covariance pairs when r = 6 in ( [T5] ), 
where Oab, acd on every panel represents the two-dimensional posterior density of (aab, o'cd)- MCMC, 
Markov chain Monte Carlo; CMC, consensus Monte Carlo; SDP, semiparametric density product; VB, 
variational Bayes; WASP, Wasserstein posterior. 


where oc has prior mean Ooc/ba- The hierarchical model for probabilistic parafac implies that 


pr(xii = Cl,... ,Xij = Cj,... 


Xip — Cp J — 7t, 


Cl ,...,Cp 


j=l 


(17) 


h=l 


.(j) 


The XijS are sampled independently given the latent class zi and probability vectors (h = 1,..., oo). 
The latent class for every sample is generated using the stick breaking representation of Dirichlet processes. 
The Gibbs sampling algorithm developed in Dunson and Xing| ( |2009| is very slow even for moderate sample 
sizes. This example demonstrates that WASP can easily scale existing sampling algorithms to massive data, 
even when efficient VB alternatives are unavailable. 

We followed the simulation setup in Dunson and Xing ( |2009 1, except with a much larger sample size. 
We fixed the sample size, number of dimensions, and number of categories in each dimension at n = 
10^, p = 20, and dj = 2 (j = 1,... ,p), respectively. These choices of n, p, and djs ensured that 
computations for sampling from the full data posterior were tractable. Data were simulated as a mixture 
of two populations such that any sample belonged to the two populations with equal probability. The two 
categories in every dimension excluding 2,4,12, and 14 were simulated from a discrete uniform in both 
populations. The dependence across dimensions 2,4,12, and 14 was induced as follows. The probabilities 
7t2,714,7ti2j and 7114 Were set to (0.20, 0.80), (0.25, 0.75), (0.80, 0.20), and (0.75, 0.25) in the first population 
and to (0.80, 0.20), (0.75, 0.25), (0.20, 0.80), and (0.25, 0.75) in the second population. The simulation setup 
was replicated 10 times. 

We used CMC, SDP, and WASP to approximate the full data posterior distributions for pr(xi = 1], 
where I G (2,4,12,14}. Two values of k G {5,10} were used for CMC, SDP, and WASP. The full data were 
randomly partitioned into k subsets and subset posterior samples for WASP were obtained after modifying 
the Gibbs sampling algorithm in Dunson and Xing (20091 using (|^. Examples for the application of CMC 
and SDP were unavailable for Dirchlet process mixtures, and it was unclear how to raise the prior density 
to the power 1/k when the prior distribution has an atomic form similar to that in ( [T^ ; therefore, we did 
not raise the prior to a power of 1/k for sampling from the subset posterior distributions in CMC and 
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Figure 4: Kernel density estimates of the marginal posterior densities for dimensions 2, 4, 12, and 14. 
MCMC, Gibbs sampling algorithm of Dunson and Xing (20091; CMC, consensus Monte Carlo; SDP, semi- 
parametric density product; VB, variational Bayes; WASP, Wasserstein posterior 


SDP The sampling algorithm for WASP based on stochastic approximation is summarized in Section |23] of 
Supplementary Material. Subset posterior samples for pr(x 2 = l),pr(x 4 = l),pr(xi 2 = l),andpr(xi 4 = 1) 
were combined to obtain their approximate posterior disttibutions using CMC, SDP, and WASP. 

The accuracy ( [TT] ) of CMC and SDP in approximating the full data marginal posterior distribution de¬ 
pended on k, with WASP outperforming CMC and SDP when k = 5 (Table |^. The approximate and full 
data posterior distributions were centered at the same value across all dimensions and replications, but the 
posterior densities for CMC and SDP were highly concentrated compared to the full data posterior density 
when k = 5 (Figure |^. The accuracy of WASP remained stable with varying k, providing an empirical 
verification of Theorems |3.1| and |3.2| in cases where our theory is not applicable. The time spent in combin¬ 
ing subset posterior samples was negligible compared to the time spent in sampling; therefore, WASP could 
be used for data with much larger sample size by choosing k large enough such that sampling was efficient 
across all the data subsets. 

Table 6: Accuracies of the approximate marginal posterior distributions for dimensions 2, 4, 12, and 14 in 
( [T^ . The accuracies are averaged over 10 simulation replications. Monte Carlo errors are in parenthesis. 
CMC, consensus Monte Carlo; SDP, semiparametric density product; WASP, Wasserstein posterior 





k = 5 



k= 10 




CMC 

SDP 

WASP 

CMC 

SDP 

WASP 

pr(x2 

= 1) 

0.63 (0.02) 

0.62 (0.02) 

0.97 (0.01) 

0.96 (0.02) 

0.95 (0.01) 

0.97 (0.01) 

pr(x4 

= 1) 

0.63 (0.02) 

0.62 (0.02) 

0.97 (0.01) 

0.96 (0.01) 

0.95 (0.02) 

0.97 (0.01) 

pr(xi2 

= 1) 

0.62 (0.02) 

0.62 (0.02) 

0.97 (0.01) 

0.95 (0.01) 

0.96 (0.02) 

0.97 (0.01) 

pr(xi4 

= 1) 

0.64(0.01) 

0.63 (0.01) 

0.97 (0.01) 

0.96 (0.02) 

0.95 (0.02) 

0.97 (0.01) 


4.5 Real data analysis: MovieLens ratings data 

We used MovieLens data to illustrate the application of WASP to large-scale ratings data. MovieLens data 
are one of the largest publicly available ratings data with about 10 million ratings from about 72 thousand 
users of the MovieLens recommender system. Each observation in the database consists of a user, movie, 
rating of the movie from 0.5 to 5 in increments of 0.5, and the time of rating. Every movie is also classified 
info af leasf one of fhe 19 genres. We fif a linear mixed effecls model ( [T5] ) using movie- and user-specific 
informafion as predicfors and fhe rafings as responses. 
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We generated three new predictors for accurate modeling of ratings following Perry (20161. First, movie 
genres were grouped into movie categories to reduce the number of genres from 19 to four: Action cate¬ 
gory included Action, Adventure, Fantasy, Horror, Sci-Fi, and Thriller genres; Children category included 
Animation and Children genres; Comedy category included Comedy genre; and Drama category included 
Crime, Documentary, Drama, Film-Noir, Musical, Mystery, Romance, War, and Western genres. If a movie 
belonged to multiple genres, then movie category scores were fractions proportional to the number of genres 
in the respective categories. Second, popularity predictor was defined as logit{(l -|- 0.5)/(n -|- 1.0)}, where 
I and n respectively were the number of users who liked and rated the movie in 30 most recent observations 
for the movie and logit (x) = log Third, previous predictor was defined fo be 1 if fhe user liked fhe 
previous movie and 0 ofherwise. We used Action, Children — Action, Comedy — Action, Drama — Action, 
popularity, and previous as fhe fixed and random effecfs in ( [T5] ). 

Following fhe sefup in Secfion 4.3 we compared fhe performance of WASP wifh CMC, SDP, and VB 
using fhe full dafa posferior disfribufion as fhe benchmark. Sampling using fhe HMC algorifhm in Sian was 
prohibitively slow for fhe full dafa posferior disfribufion, so we firsl randomly selecfed 5000 users and fhen 
randomly selecfed 20 ralings for every user. This resulfed in a dafa sel wifh 100,000 ratings. We randomly 
splif fhe users info 10 framing dafa sels such lhal ratings for any user belonged fo fhe same fraining dafa 
sef. To compufe fhe approximafe posleriors using CMC, SDP, and WASP, we sel k = 10 and randomly 
partitioned fhe users info k subsels such lhal each subsel conlained all fhe ralings for a user. This sefup was 
replicaled for every fraining dafa. 

WASP performed beffer lhan ils compefilors in approximating fhe full dafa posterior dislribufions for 
variances and covariances of fhe random effecfs. Similar fo fhe simulation resulls in Secfion |4.3| CMC, 
SDP, VB, and WASP were significanlly fasler lhan fhe full dafa posterior disfribufion, wifh VB being fhe 
faslesf. CMC, SDP, and WASP showed excellenl performed in approximaling fhe full dafa posferior dislri- 
bulions for fhe fixed effecfs. WASP oufperformed ils compefilors in approximating fhe full dafa posferior 
dislribufions for variances, covariances, and pairs of covariances of fhe random effecfs (Tables [7| and 
1^. VB significanlly under-performed in fhe eslimalion of fhe posferior disfribufion for fhe fixed effecfs and 
covariance malrix of fhe random effecfs. The accuracy of marginals in CMC and SDP depended on fhe mag- 
nilude of covariances, wifh bolh melhods showing excellenl accuracy for covariances wifh low magnilude. 
The accuracies of fhe Iwo-dimensional join! dislribufions in CMC and SDP were poor because fhe full dafa 
posteriors concenlraled af differenl locations (Figure |^. Excepl for fhe poor performance of CMC, SDP, 
and VB in approximating fhe posterior disfribufion of variances and covariances of fhe random effecfs, our 
resulls for MovieLens dafa agreed wifh fhe simulafion resulls. We concluded lhal WASP performed heller 
lhan ils compefilors in fhe analysis of MovieLens dafa. 


Table 7: Accuracies of fhe approximafe posleriors for variances in ( [T5] ). The accuracies are averaged over 
10 replications. Monfe Carlo errors are in parenlhesis. CMC, consensus Monle Carlo; SDP, semiparamefric 
densify producl; WASP, Wassersfein posferior 



0-2 

Action 

0-2 

Children — Action 

0-2 

'^Comedy — Action 

0-2 

Drama — Action 

fT^ 

Popularity 

0-2 

'^Previous 

VB 

0.00 (0.00) 

0.00 (0.00) 

0.00 (0.00) 

0.00 (0.00) 

0.00 (0.00) 

0.00 (0.00) 

CMC 

0.28 (0.13) 

0.01 (0.01) 

0.01 (0.01) 

0.14 (0.09) 

0.74 (0.10) 

0.22 (0.10) 

SDP 

0.05 (0.03) 

0.00 (0.00) 

0.00 (0.00) 

0.01 (0.01) 

0.35 (0.10) 

0.03 (0.03) 

WASP 

0.92 (0.04) 

0.93 (0.02) 

0.87 (0.06) 

0.85 (0.08) 

0.92 (0.03) 

0.93 (0.05) 


5 Discussion 

We have presented WASP as an approach for compulalionally efficienl approximation of fhe posterior dislri- 
bulions of paramelers and Iheir funclions when fhe sample size is large. WASP allows exfensions of existing 


16 









Table 8: Accuracies of the approximate posteriors for covariances in ( [T5] ). The accuracies are averaged 
over 10 replications. Monte Carlo errors are in parenthesis. The subscripts 1,..., 6 are used for predictors 
Action, Children — Action, Comedy — Action, Drama — Action, popularity, and previous. CMC, consensus 
Monte Carlo; SDP, semiparametric density product; WASP, Wasserstein posterior 



0 ’l2 

ffis 

0'l4 

ffl5 

Cie 

0’23 

^■24 

0'25 

VB 

0.00 (0.00) 

0.00 (0.00) 

0.00 (0.00) 

0.00 (0.00) 

0.00 (0.00) 

0.00 (0.00) 

0.00 (0.00) 

0.00 (0.00) 

CMC 

0.06 (0.03) 

0.16(0.04) 

0.18(0.04) 

0.83 (0.07) 

0.33 (0.13) 

0.01 (0.01) 

0.07 (0.02) 

0.80 (0.04) 

SDP 

0.01 (0.01) 

0.08 (0.03) 

0.07 (0.02) 

0.75 (0.06) 

0.14 (0.09) 

0.00 (0.00) 

0.02 (0.01) 

0.73 (0.08) 

WASP 

0.95 (0.02) 

0.91 (0.04) 

0.91 (0.05) 

0.94 (0.03) 

0.90 (0.07) 

0.89 (0.07) 

0.85 (0.08) 

0.93 (0.03) 


0’26 

0’34 

^35 


Crs 

0'46 

tTse 


VB 

0.01 (0.00) 

0.01 (0.00) 

0.01 (0.00) 

0.01 (0.00) 

0.00 (0.00) 

0.00 (0.00) 

0.00 (0.00) 


CMC 

0.66 (0.09) 

0.65 (0.07) 

0.76 (0.08) 

0.71 (0.05) 

0.82 (0.04) 

0.61 (0.11) 

0.55 (0.09) 


SDP 

0.59 (0.11) 

0.62 (0.06) 

0.64 (0.09) 

0.66 (0.08) 

0.66 (0.09) 

0.56 (0.14) 

0.55 (0.13) 


WASP 

0.91 (0.05) 

0.94 (0.05) 

0.93 (0.03) 

0.91 (0.04) 

0.93 (0.04) 

0.93 (0.04) 

0.94 (0.04) 



Table 9: Accuracies of the approximate two-dimensional joint posteriors for the covariances of random 
effects. The accuracies are averaged over 10 replications. Monte Carlo errors are in parenthesis. The 
subscripts 1,..., 6 are used for predictors Action, Children — Action, Comedy — Action, Drama — Action, 
popularity, and previous. CMC, consensus Monte Carlo; SDP, semiparametric density product; WASP, 
Wasserstein posterior 



(o'i2, nia) 

(o'i2, ffir) 

(0'l2, ITls) 

(o'i2, ffie) 

VB 

0.18(0.04) 

0.22 (0.07) 

0.31 (0.03) 

0.31 (0.02) 

CMC 

0.05 (0.02) 

0.04 (0.02) 

0.06 (0.03) 

0.05 (0.02) 

SDP 

0.05 (0.02) 

0.04 (0.02) 

0.06 (0.03) 

0.05 (0.02) 

WASP 

0.88 (0.03) 

0.88 (0.03) 

0.88 (0.02) 

0.86 (0.06) 


samplers to massive data with minimal modifications and is easily implemented using probabilistic pro¬ 
gramming languages, such as Stan. Theoretically, we showed that the rate of convergence of WASP to the 
delta measure centered at the true parameter value in W 2 distance matches the optimal parametric rate up to 
a logarithmic factor if the number of subsets increases slowly with the size of the full dataset. Empirically, 
we demonstrated that results from WASP and MCMC agree closely in several widely different examples, 
while WASP enables massive speed-ups in computational time. 

We plan to explore several extensions of WASP in the future. First, it is unclear how to extend stochastic 
approximation to cases where the likelihood is unavailable in a product form. This extension in crucial for 
proper uncertainty quantification outside of settings in which the observations are conditionally independent 
given latent variables. Second, it is unclear how to optimally choose k in practice; larger k improves 
computational time when abundant processors are available but choosing k too large may lead to increasing 
statistical errors (refer to Theorem 3.11. Our numerical experiments show that the accuracy of WASP is 
robust to the choice of k if all the subset sizes are moderately large relative to the number of parameters. In 
addition, it is of interest to study more deeply the impact of the partitioning schemes and attempt to develop 
approaches that deal with not only large sample sizes but also high-dimensional data. A possibility in this 
regard is to combine WASP with approximate MCMC ([Johndrow et al.[[2015p. 
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Figure 5: Kernel density estimates of the posterior densities of four covariance pairs, where ffab, Ccd on 
every panel represents the two-dimensional posterior density of (o’^t,) o'cn)- MCMC, Markov chain Monte 
Carlo; CMC, consensus Monte Carlo; SDR, semiparametric density product; VB, variational Bayes; WASP, 
Wasserstein posterior. 


Appendix 


A Assumptions 

Suppose h,(pi,p 2 ) = [J(y^ — y^)^d'v] is the Hellinger distance between two generic densities 
Pi,P 2 - We define the following pseudo Hellinger distance on each subset of data Y[j], with sample size m 
and j = 1,..., k. 

Definition A.l (Pseudo Hellinger distance) The pseudo Hellinger distance between and Pg^ 

in 1 ^ 0 , 3 ,i : 6 £ 0} is defined as 


I I ®2)) 


(18) 


i=l 


This definition is based on Equation (3.1) in Ghosal and van Der Vaart (20071 and it is straightforward to 
check that for every j = 1,..., k, ({®[^iP 0 ,j,i : 9 £ ©}, h.Tn,j) is a metric space. 

For a m-dimensional random vector Z = (Zi,..., Zm)^, denote its Lq norm as 


|Z|q = 


1 1/q 


m 




jq’ 


i=l 


and we also use ||Z|| to represent IZb- 

Next we define the generalized bracketing entropy for inid data. If the data are indeed iid then the 
generalized bracketing entropy coincides with the usual bracketing entropy. 

Definition A.2 (Generalized bracketing entropy) Let Z be a fixed subset of&. For a fixed j £ {1,..., k}, 
let 


Tj(Z) = {pj(y|0] = (pji(yi|0),...,Pjm(yml0))'^ -V = (yi,... £ ®lliyji,0 £ Z} 

be the class of m-dimensional functions indexed by 0. For a given 5 > 0, let 

® ^s] • 1-s (y ] = (^-sl (y l)) • • • ) l-STa(y m)) ) IJ-s (y ] = (y l) j • • • ^ trsmty m)) ) 
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be the generalized bracketing set of7^['E) with cardinality N, such that for any Pj(y|0) G J*] there 
exists a pair of functions [Ig, Ug] G 23 (5, J*] such that 


l-si(yi) ^Pjt(yi) for ally G and alii = I,... ,vn 

and ||\/ul- V^ll ^ 

The -bracketing number of N[] (6, ^Pj (E), is defined as the smallest cardinality of 

the generalized bracketing set 23 (6, (El)). The h-^^j-bracketing entropy of IPj (El) is defined as 

H[] =log(l + N[] 


B Proofs of Theorems 

/ \-l/(2a) 

Let em, = ( lo^Ta J ■ notation, in all the following proofs, we will sometimes write 

p(yji I 0) =Pji(yjt I 0]- 


B.l Proof of Theorem 13.II 

Due to the compactness of © in (Al), we assume that p(0, 0o) ^ Mq for a large finite constant Mq. We 
start with a decomposition of the W 2 distance from the jth subset posterior n^^(- | Y[j]) to the delta measure 
at the ttue parameter 0o: 


EpeoW|(n];(- I Y[p),5eJ-)) =EP3^ p2(0,0o)^^;(d0 | Y[p) 


0 


p2(0,0o)nT;(d0|Y[p) + Ep3^ 

J{0:p(0,0o)^Ciem} 

+ M^Ep n];, (p(0, 0o) > Clen^ | Y[j]) . 


{ 0 :p( 0 , 0 o)>ciem} 


p2(0,0o)nx,(d0|Y[p: 


(19) 


We will choose the constant Ci as Ci = ( 


l/[2a] 


, where gi, Cp, qi, ti are the constants in (Al), (A2), 
and Lemma [T75] and Lemma 1.6 in Supplementary Material. 


The following proofs are similar to the proofs of Theorem 1, 4, and 10 in Ghosal and van Der Vaart 


(20071. The main difference is that our likelihood has been raised to the power y. Using condition (A2), we 
can further replace the p metric by the pseudo Hellinger distance: 

(0 G 0 : p(0, 0o] > Ciem I Yfp) 

^ (0 G 0 : H^j(P 0 ,j,P 0 o,j) > | Y[p) 


n™i 

r p(Ynl0) 

Lp(Vjx|0o)J 

^n(d0) 

{0e0:h^P0,0o)>y%if efc} 

r p(Ynl0) 
Lp(Yji|0o)J 

^n(d0] 


( 20 ) 


For the denominator in ( |20l ), by Condition (A4) and Lemma [L6l for m sufficiently large, with probability at 

,2a) 


least 1 — exp(—r 2 Tae:J^ 


Ta 



p(Yjt|0F 

p(Y3t|0o)Y 


n(d0) > exp(—TiTie^). 


( 21 ) 
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For the numerator in ([20|), by Condition (A3) and Lemma 1.5 we set 6 = Y^ 2 Tig 2 /qie^ and obtain 


that with probability at least 1 — 4 exp ( — j ^ 1 — 4 exp ( — ), 


.?Il3A-n c-2a 


sup 


m 

rr 

\ p(Yjiie] 1 

11 

Lp(Yii| 0 o]J 


^ exp (- 2 Tig 2 me^) ^ exp (- 2 rine^) ( 22 ) 


Therefore, based on (| 2 T), and (| 22 l), we obtain that with probability at least 1 — 


exp(-r 2 me^), 


4 exp 


Let Ae^ be the event {0 G © : IT ^0 G 0 : p(0, 0o) > Ciera Y[j] j ^ exp (—Tine^) }. Then we can 
bound the second term in ([T^ as 


(^0 G 0 : p(0, 0o) > CiEra Y[J]^ ^ exp (-2rine^ + rine^) ^ exp (-tine^) . 


tPeo'^m (p(0,6o) > ClOra 

[l(Ae,,)n^;(p( 0 , 0 o] 


Y[j] 


> CiEt 


Y[i] 




(p( 0 , 0 o)>ci( 


Y[j] 


^ exp (-rine^) + P 0 ^^(A^^) • 1 
^ exp (-nne?^) + ■ 

^ 6 exp (-C 2 me^) , 


^ exp (—Tine^) +4exp ( —) + exp(—r 2 me^‘^’' 


for C 2 = min(ri, r 2 , 2 riq 2 /qi], as clearly the second term is dominating the other two given m < n. 
Therefore, for ( [T^ , since era = (m./ log^ as m —)■ oo, an explicit bound will be 

EpeoW|(n];(- |Y[p),5eJ-)) ^ ^ + 6 Mg exp (-C 2 log^ m) 


2 log^/“ m 1 _ 

^ ^ Cl 

^ J- _l/rv l_i_4^ 


TU 

logA« m 


m^/“ m^/°^ 

as m becomes sufficiently large, where the constant Ci depends on a,Ci,C 2 , which further de¬ 
pends on gi, g 2 , qi, q 2 , Pi,P 2 ) Cl- Since qi,q 2 in Lemma 1.5 and ri,r 2 in Lemma 1.6 depend on 
gi, g2, Di, D2, K, Cji, it follows that Ci depends on gi, g2, Cl, Di, D2, k, c „. □ 

B.2 Proof of Theorem 13.21 

Based on Theorem|3.1|and Lemma 1.7 in Supplementary Material, we have that for any C > 0, 


^W2(n^MY(’')),S8 ,(-)) >y'^ 


C log^/'^ m 


l/a 


(n) 

I k 


C log^^”^ m 


(L 1 


C m *^^0 


n 2 


i^W2(n^;(.|Y[p),50j.)) 


i=i 
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(ii) Jl 

log^/^^m _ Cl 
^ Cklog^/^^m ^ Tai/« C’ 

where we have used Markov’s inequality in (i) and the relation between li and I 2 norms in (ii). This implies 
thatW2(TT;(.| Y(-)),5eJ0) =Op □ 


C Univariate density estimation 

Let Xi,... ,Xn be n copies of a scalar random variable X that follows probability distribution Pq with 
density po- The full data are randomly split into k subsets and Xji,..., Xjm, represent the data on subset 
j (j = 1,..., k). The hierarchical model for density estimation using the stick-breaking representation of 
Dirichlet process mixtures is 

00 

Xji I Zji,{iJ,h}h=i,{o'H}H=iZji~^v/H6H, "‘'H = Vh (1 - Vi), Vh | ~ Beta(l, «), 

h=l t<h 

a ~ Gamma(a(x, hoe), pn I cTh ~ ^h)) ~ Inverse-Gamma(ao-, bo-], (23) 

where Urr > 2 and Beta, Gamma, and Inverse-Gamma random variables have means and 

L~t“OC O (X ^(T J- 

and variances y.;-——v, and -yrr. If I* is the maximum number of atoms in the stick- 

breaking representation, then the prior density n is in the form a discrete mixture. We cannot use existing 
sampling algorithms directly if n is raised to a power of 1 /k, so it is unclear how to sample from the subset 
posterior density of competing approaches in Section [2^ 

We show that it is still possible to sample from the subset posterior density in Q using data augmenta¬ 
tion. Let Lj be the likelihood given Xj i,..., Xj m, and latent variables zj i,..., Zj m, in ( [^ , then 




.“2i (tm 

2 e ’, 


(24) 


h=l 


where l(zjt = h) is 1 if Zjt = h and 0 otherwise and (Jhj = l(zjt = h.). For stochastic approxima¬ 
tion, we raise Lj in ([24l) to the power y and obtain 


lYcr it* r 2 it* r it* i 1 Tfo 2-1 — 2^21 IIi=i M-n) YttH, 

1-]^({Ph}h=i,{o-h}h=ii{^h}h=i) = li(27raH) ^ e '"h v/^ k 


h=l 


Standard arguments imply that the analytic form of full conditional densities of parameters are 


Ph I rest oc e 


2a4 Ytthj+1 J 


vny _ _ Viifv _ii,,F 1 —-bh , _i>CT 

cT^^lrestoco-^^^e Li.ii(z,-h)(x„ M ^ 2 - 2 ^ e 

Vh I rest oc Lz,=h) (1 _ 

t* 

a I rest oc ]^(1 —Vd)‘ 


\ a—1 


h.=l 


(25) 


(26) 
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for h = 1 


...,r.Let 


^ yLi=iUzji = H)xji < 


Q-jh. — 


ytJHj +1 

ytjHj +1 


m 2 

+ Oa, bjH = ^ X l(Zji = H) (Xji - |Xh)^ + ^ + be 


2 ’ J ' ^ 2 rii-/ 2 

i-1 

for h, = 1,..., I*, then all full conditional densities are tractable in terms of standard distributions: 
jj-j H I rest ~ N (mj h,, Vj vi ], cr?|^ | rest ~ Inverse-Gamma (Oj h , bj n) j 

m m 

Vjh, I rest ~ Beta(l + y L l(zji = h],a + y^l(zji >li)), 

t* 


1 = 1 


1=1 


ajH I rest ~ Gamma(acc + r,b(x - ^ log(l - VjH))- 


h=l 

Finally, posterior distribution of the latent variables is 

PjH = “ 

h=l 


Zji I rest ~ Z Piti^H) Pj 






(27) 

(28) 


(29) 


(30) 


where vjh = Vjh Ot<H(3 ^ ^jt) [N’(Ta, v) is the Gaussian density with mean m and variance v. 


D Linear program 


minimize 

a,Ti,...,T]e 
subject to 


k 

y~ trace(Tj^Dj) 
j=i 


0 ^ Gi ^ 1, i = l,...,g, 

rr = l,...,g, v = l,...,sj, j = l,...,k, 

1^ a = 1, 

Tjlsj=a, j = 

T 71 s = ^, j = l,...,k. 

Sj 


(31) 


This linear program can be solved using a variety of linear programming solvers in Mat lab or R. More 
specialized algorithms are developed in Cuturi and Doucet (|2014|),[Carlier et al.|(|2015|), and|Srivastava et ah 


(20151. We provide a simple implementation using Matlab in Sectionj^of Supplementary Material. 
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Supplementary Materials for Scalable Bayes via Barycenter in 

Wasserstein Space 


1 Technical Lemmas 


To show Theorem 3.1, we first introduce in Lemma [T31 a generalized version of the concentration in¬ 


equality in Theorem 1 of Wong and Shen (19951. The proof parallels the original proof in Wong and 


Shen (jl995 1, with several adaptations for the inid setup. Let Zji(0) = log[p(Yji|0)/p(Yjt|0o)], and 


Zji(0) = max(Zjt(0), — t) be the lower truncated version of Zji(0) for some constant t > 0 to be chosen 
later. Let Zj (0) = (Zji(0],... ,Zj^(0))T and Zj(0] = (Zji(0),..., Zj^(0))T. 


Lemma 1.1 Let Cit = 2e — e Then for any 0 G 0, 


1 

-}^EZjt(0) ^-(l-ci,)h2^^(0,0o). 


(32) 


Proof of Lemma im 

The proof is a simple adaptation of Lemma 2 and Lemma 4 in [Wong and Shen] ( |1995| ). First note that for 
every observation Yjt (which takes value rjji), we can define fhe event Aji = {yjt : p('yji)|0)/p(yii|0o) < 
e“^}. Their Lemma 2 implies that P(Ajt) ^ (1 — e“'^/^]“^h.^(pji(-|0),pji(-|0o]), which further implies 
by simple averaging over i = 1,..., m that 

1 

- Y. ^ (0,0o). (33) 

i=l 

Following the same derivation of their Lemma 4, we have for every individual Zjt, 

EZji(0) ^ -H2(p3,(.|0),pj,(.|0o]) + 


Then a simple averaging over r = 1,..., m together with (331 gives (32). □ 


Lemma 1.2 Let C 2 t = — 1 — t/2)/(1 — e For any t > 0, integer f ^ 2 and any 0 G 0 that 

satisfies h-Tn,] (0, 0o) ^ T 


m ^ 

i=l 


Zji(0) 


2V2C2t:T 


f! f 1 


2 VV2C2t1’ 


1-2 


Proof of Lemma [L2t 


Lemma 5 in Wong and Shen (19951 is stated for every single observation Yjt, so it still holds for individual 
Zji(0). For all 1 = 1,..., m, 


So 


exp 


Z3t(0)/2 - 1 - Zji(0)/2 < C2xh^(pji(-|0),P3i(-|0o)), 


where C 2 t is as defined in the statement of the lemma. Averaging over i = 1,..., m yields 

m 


1 


m 




exp 


Zji(0)/2 -1- Zji(0)/2 


i=l 


^ C2TfL?n,j(0,0o) ^ C2T■r^ 
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where the second inequality is from the condition hm.) ( 6 , So) ^ t"- Using — 1 — x ^ for all x > 0, 
we have 


1 


m 




i=l 


^ 2HlC2rr. 


Rearranging the terms and the conclusion follows. □ 

Lemma 1.3 Let j G {1,..., k} be fixed. Suppose E is a subset ofO. Let Zj (E) = |Zj (0), 0 G !e| . For 
u > 0, 


any 


H 


n (u,£j(Z),|| • ll) ^ . 


Proof of Lemma [131 

The proof follows the argument in the proof of Lemma 3 in Wong and Shen (19951. We can derive that for 
each r = 1,..., m, for any 0i, 02 G Z, 


-Peo 


1 2 


Zjt(0i)-Zjt(02)J ^4e"h^^^(0i,02). 

Then averaging over i= 1,..., m gives the relation between two norms 

Zj(0i)-Zj(02)|| ^2e"/2^^j(0i,02), 

which further implies the relation between the bracketing entropies. □ 

Lemma 1.4 ( van der Geer and Lederet\ {2013) Theorem 8) Let ] G {1,..., k} be fixed. Suppose a class of 
functions Tj = {f(y) = (fi(yi), • • •, fm(ym))^,R = (Ui, • • ■ ,'yTa) G (8)DLi Dji} satisfies 

(i) supfg5-. ||f|| ^ 1; 

(ii) For any integer I ^ 2, supf^gr. |f|^ ^ /2, for some constant M. > 0; 

Then for any t > 0, 


( .j III. 

Y. “ tp fi(Yji) 


^ min Rs + — + 24\/^) ^ 2e~^ 

SeN VTTL 


where 


Rs =2-SVni+14\/6^2-YH[] (2-s,Tj, 


+ 


36MHu (l,3-j,||-||) 


s=0 


m 


Proof of Lemma lL4t 

The theorem we present here is slightly different from the original Theorem 8 in [van der Geer and Lederer 


(20131 in that we have used the generalized bracketing entropy in Definition A.2 in the manuscript. Although 
the original version is presented with the usual L 2 -bracketing entropy for univariate functions, the whole 
“chaining along a tree” argument will still go through, if we replace the || • || and | • |q norms and bracketing 
entropies in their proofs by our generalized versions to multivariate functions as in Definition A.2. □ 


Lemma 1.5 { Generalization of Wong and Shen ( 1995 1 Theorem 1) Assume (A3) holds. Then for any 5 > 0, 
there exist positive constants qi, q 2 that depend on Di, D 2 , such that for all subsets Yjp with j = 1 ,..., k 
and all sufficiently large m, 


)(n-) 


sup 


n 


p(Y)i I e) 


H,,j(e,0o)>6fL|p(Yji I 0o) 


^ exp(—qimS^) ^ 4exp(—q 2 m 6 ^ 


(34) 
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Proof of Lemma [131 

We consider the class of functions 


Zj(r) = 


f Zj(0) 


: 0 G 0 satisfies h,Tn,i ( 0 , 0 o) ^ t , 


I 2^/2^r 

for a fixed r > 0. This class is a rescaled version of 

Zj ({0 G 0 : hraj (0, 0o) ^ f}) = |Zj (0] : 0 G 0 satisfies (0,0o] ^ r| 


as in Lemma 


1.3 


By Lemma|l.2[ 2.j(r) satisfies Conditions (i) and (ii) in Lemma 


1.4 


wifh fhe consfanf 


M = l/(-v/2c^r). Therefore fhe concenfrafion inequalify in Lemma ] 1.4| holds for Zj (r). 

We firsf simplify fhe ferm Rs in fhe inequalify. Rs involves fhe L 2 -brackefing enfropy of fhe class Zj (r). 
is nonincreasing in u, we have 


Since H[] [u, Zj (r), 
s 


^2-^Jhu (2-^^j(r),||•||) ^2^ 

s=0 ^ s=0 


S .2-^ 

2 -(s + l) 


H[] (u,2;j(r], 


du 


= 2 


= 2 


2-(S+l) 

■1 

2-(S+i) 


H[] u,2,j(r), 


du 


(^2\/2c2Tru, ({0 G 0 : h,raj(0,0o] ^ t}), 


Hr 


1 


(ii) 


\/2C2tT' 
1 


v/2c2^r 


■2V2C2tT 

2-S^2C2tT 
'2v'2C2t1' 

2^5 V2C2t1' 
%/2c^4e^r 


du 

du 


H[] (^u, Zj({0 G 0 : h,TTT,j(0,0o] ^ r}), || 

(u/(2e^/2),(Pj({0 G0:h^j(0,0o) ^r}),li^j)du 


2-(S+i)y'2c2Te-'"T 


(u,Tj({0 G 0 : h,n^j(0,0o) ^ r}), h^aj) du, (35) 


where (i) follows fro m fh e scaling relation befween Zj(r) and 2.j({0 G 0 : h,m,j(0,0o) ^ t})! and (ii) 


follows from Lemma 


1.3 


We choose infeger S ^ 1 such fhaf 2-(s+2) ^ 72c2Te-"'r/2i2 ^ If 

is possible fo do so because we only need fo consider r ^ \/2 (since fhe hraj disfance is upper bounded by 
\/2), 02 x 6 “'^ ^1/2 for all t ^ 0, and if is guaranfeed fhaf A/ 2 c 2 ^^e^r/ 2 ^^ ^ y/2/2^‘^ < 1/4. 

Now we can apply (A3) fo ([35]) and obfain fhaf uniformly over all j = 1,..., k, 




s=0 

^ . - • D2VtTL 




\/2c2Te 

2 A/'k(u, r)du 

(V2c2Te-^r] /2i2 


V2c2Te-^ \ 2D2\/2c2Te-^ 

T 1 =-ztt; -Vrar. 


Di 


D? 


(36) 


Furfhermore, since (A3) says V(u, r) is nonincreasing in u, we can also derive fhaf 


H 
^ H 


U (l,£j(r),|| • ll) =H[] (2V^r,Zj({0G0:h^j(0,0o) 

[] (^y/2c2Te-'^r, Tj ({0 G 0 : h.TTT,j (0, 0o) ^ r}),hTaj) ^ (^V2c2Te-^r,r^ 
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1 


V2c2Te-^r - (x/2c2^e^r)^2^2 


\/2c2Te '^T 

(\/2c2Te-'^T)^/2i2 
T ^^2 


n 2 


-^^(u, r)du 


2D|c2Te '^mr^ ^ 8D|c2Te ^ttlt 

2 < 54 


Df (1 - V2c2Te-^r/2i2) 
s 

By our choice 


( 37 ) 


■^1 


where in the last step, we used the fact that \/ 2 c 2 re~'^r/ 2 ^^ < 1/2. 

By our choice 2-(S+2) ^ ^2c2Te-'^r/2i2, it follows from @ and ([^ that 


minRs ^ 
seN 


V2C2re-^ 


210 




V2 


Df 


36 


'— . ^ /7^ 2D2\/2c2Te /— , 

mr + 14v 6 • ^/mr H— . 

V2c2T'nir 


8D|c2Te "^mr^ 


Df 


_ +56^3^ +144 V2^e-"/2 


\/ C2i 


(38) 


In the inequality of Lemma 1.4 we let t = Ca^mT^ with integer £ ^ 1 and constant 03^ to be chosen later, 

I 2 

then with probability at least 1 — , the empirical process 


^ m 

~7^ X - Ep0^fi(Yji 

fezd’-) ^ i=i 


will not exceed the following upper bound 

_L _1_ 133,/oE 

210 


min Rq + 
seN 




^ + 56^/3^ + 144^/23e-"/2 

Df 


'D| 


V C2^:e-'^^/mr + + 2A^/Qc^^r 

V2c2T'm.r 


36 

: C4T:V'nir + , 

V2c2TTnr 


where C 4 t: = ^ + 56\/3^ + 144\/2^e \/c2r£~'^4- + 24-y/6c3^. As a result, the empirical 

" U V'"^2t 

~ 2 
process on the class Zj ({0 G 0 : hiTij (0, 0o) ^ r}) satisfies that with probability at least 1 — 2e^‘^3Tmr ^ 


sup 

0e0:hn^j (0,eo)^r 


/m f 


X [Zji(0]-Epg^Zji(0)] ^ 2V^r- (^CirV^r + 


72 

^ 2V2c2tC4t • \/mT^ H— 

Vm 


(39) 


On the set {0 G 0 : r ^ hm.) (0, 0o] ^ 2 t|, we have h^- (0, 0o) ^ r^. By Lemma 

1 

— y* EZji(0) ^ -(1 -Cit-)t“. 


1.1 


it follows that 


m ■ 


i=l 


This together with ( [39] ) implies that with probability at least 1 — 2e 

m 


-caTmr* 


1 ~ 72 

sup — Y Zii(0) ^ - (l - CiT - 8\/2 c2tC4t) H-• 

0e0:r^h„,j (0,0 o)^2t TTl TTL 
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Now we choose t such that e = 1/32. Then Cix < 0.07, 29 < C 2 t < 30. Set Csx = 1/2^°. By (A3), 

D2 ^ D?/2^2_ 3o 


8\/2c2'tC4T ^ 8\/60 • 


V2 , 56\/3 , 144^2 
^ + 16.224 


, 30 


36 


210 ^. 230 


+ 24i/^!>< 0.377, 


which implies that 1 — Cit — 8-v/2c2tC4t > 0.55. Thus we have proved that for any 0 < r ^ \/2, uniformly 
for all j = 1,..., k and for all sufficiently large m. 


(n) 


sup 


{0e0:rs;H,nj(0,0o)sS2r} TU ^ 


1 ^ 72 

Zjt(0) ^-0.55r2 + - 

m ^— m 


^ p 


(n) 


1 


2 . 72 


sup —y Z.]i[Q) ^ —0.55r^ + 

y{0e0:r^h,nj (0,0o)sS2r} ^ TTL 




2 exp (-mrV2^°) . (40) 


Finally for a given 6, we set r = 2^6 for integers £ ^ 0, m > 72/0.05/6 = 1440/5, and let L be the smallest 
integer such that 2'^5^ > 2. It follows that 




{0e0:h„,j(O-Oo)>6}fj^ I ®o) 


n 


p(Yii I e) 


^ exp(—0.5m6^ 


^ P 






(n) 
00 

L 

Lp 


' 1 ™ 72 

sup — y_ Zji(0) ^ —0.555^ H- 

y{0e0:h,,j(O-Oo)>6}llli ^ 


i=l 


» 


1 


f—Q ° V{0e0:2«6^H^j(O-Oo)^2«+l6} TTL ^ 
y~ 2 exp (—2^^m6^/200) ^ 4 exp (—m6^/200) 


Ta 

^Zji(0) ^-0.55(2^6) 


2^72 
m 


We set qi = 0.5, q 2 = 1/2^0 and complete the proof. 


□ 


Lemma 1.6 Assume (AS) holds. Then for any 6 > 0, there exist positive constants ri,r 2 that depend on 
gi, 92 , K, Cji, such that for every subset Y[j] (j = 1,..., \i),for any t ^ e^, 




p(Yjii0r 

p(Yjt|0o)Y 


n(d0) ^ exp(—Tint) 


^ exp(—T 2 mt]. 


Proof of Lemma [L6t 

Define fhe evenf 0e^ as 


0. 



p(Yji|0o) \ 
p(Yjil0) ) 


1 ^ e 


2a 


fhen (A4) can be wriffen as n(0e^) ^ exp(—Cyttie^). For A C 0, lef ne^(A) = n(An 0e^)/n(0e^) 
be fhe prior measure FI resfricfed measure fo fhe sef 0e„^- For fhe leff-hand side of fhe conclusion, we have 


p(Yi) 


0 


n 

i=l 


p(Yji|0r 

p(Yji|0o)^ 


n(d0) < exp(—Tint] 
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(i) 

^ P 


(n) 


(ii) 

^ P 


(tv) 

^ P 


0 . 


fr p(Vji|9F 

yp(Yjt|0o)T' 


TT(d0] ^ exp(—Tiut) 


^ p';^ n(0e. 


0 


£m \ 


p(Yji|0r 
fj^p(Yji|0o]^ 


n 


r[e^(d0) ^ exp(-rint) 


(n) 


0. 


n 


p(Yjii0r 




rfe^(d0) ^ exp(—Tint + CytTLe 


2(X\ 
m J 




(n) 


vi=l 

m 


^ Tigimt-c^game?^ 


L 


i=l L 


1 p(Yji|0o)^ ^ 


0e 


, p(Yji|0o)^ 


-l\^^2a 


^ Tigirat-(C7tg2 + K )Tae 


(41) 


In the derivation above, (i) follows from making the region of the integral smaller; (ii) is from (A4); (iii) 
is an application of Jensen’s inequality and uses the fact that giym ^ n ^ g 2 yTu; (iv) follows by using 
Fubini’s theorem, the inequality x ^ (e'^^ — 1 )/k for x ^ 0, and (A4): 


m 


Le, 

i=l 


0 


, p(Yji|0o)^ ^ 

log ne„,(d0) ^ 


em 

m r 


p(Yp|0) 


0 


Z ^Peo log_ 


i=l 


p(Yji|0o) 

p(Yp|0) 


ne„,(d0) 


0 


p(Yji|0o) 


-L ^ ^ J 1 r I ‘ 


-1 


ne^(d0] ^ K ^me 


-1^^2a 

m • 


Now in ( |44] ), if we choose ri large enough and let Zi = Jq log 1~1 em (’ "'o oan use the 


Bernstein’s inequality (Corollary 2.10 in Massart (20031) to control the large deviation for the sum of Zi’s. 
By inspecting the conditions for this inequality, we can see that for any integer f ^ 2, by the inequality 
(kx) V^- ^ ~ 1> (A4), and the Jensen’s inequality. 


£e>-.. [(Zt) 


i=l 


f! 


m 

Z ( 


kEt 


(Zi)+ -1 


V. 

^ —rtne 


2k 
1 ">.v-Ta • 


i=l 


Therefore in the Corollary 2.10, we can set c = k~^, v = 2K“^me^, x = tigimt — ( 0^02 + K~^)me^. 
If we choose ri > (C 7 tg 2 + 3K~^)/gi, then since t ^ e^, we have x > 2K~^mt ^ 2K“^Tae^. Hence it 
follows that v/c = 2K“^Tae^ < x. We apply the Bernstein’s inequality to ([44]) and obtain that 




)(n-) 


exp 


' m 


i ^Peo^i 


^ Tigimt-(C7tg2 + K £me^ 


a=i 


gexp(-^)<exp(-^^)<exp(-^^) 


2(v + cx) 

We set r 2 = 1/2 and complete the proof. □ 

Lemma 1.7 Let y denote the W 2 barycenter o/N measures wi,..., wn in (P2(0)- Then for any 0o G 0, 

1 ^ 


W 2 (’V, 60 (,) ^ W2(Vj, S0„’ 


i=i 
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Proof of Lemma [IJ} 


Theorem 4.1 in 


Agueh and Carlier (20111 shows that y = TJjYi, where 

1 




i=i 


for any 0 G 0 and is the optimal transport map that pushes yi forward to Vj, i.e. "Vj = is 

the identity operator. We use the property of the p metric in (A5) and obtain that 


W2(Y,50J = 


N 


0 




j=i 


0 


N 


n 2 


N 


Xp(T^'(0),0o) 


j=i 


Yi(d0) 


1 

W 


N If 

Y_ p2(T/(0),0o)A/i(d0) + ^}^ p(T/(0],0o)p(T,^(0),0o)w(d0]. 


(42) 


j=i 


0 


For each term in the second summation, we apply the Cauchy-Schwartz inequality and have 
p (-1^1(0), 0o)p (1^(0), 0o) W(d0) 


0 


^p2 (i:M0),0o) w(d0) • Jj^p' (T|(0],0o) Vi(d0) = W 2 (Vj, 6 eo) • W 2 (v;, 50 j. 


Therefore in (42i, 


N 


W|(v,50j ^ ^^W|(wj,50j + ^^W 2 (vj, 60 j • W 2 (vt, 60 Q; 

j=i i/t 

H 2 


1 ^ 

^^W 2 (vj, 60 Q; 


i=i 


and hence the conclusion follows. □ 


2 Experiments 

2.1 Simulated data: finite mixture of Gaussians 

Consider the set of L mixture of Gaussians. If y ^ G is the ith observation (i = 1,..., n) sampled from 
a mixture of L Gaussians, then 


L 

PtUi I 0) = I (43) 

t=i 

where 7t = (tti, ..., ttl) lies in the (L — l)-simplex, \xi and Z; (1 = 1,..., L) are the mean and covari¬ 
ance parameters of a p-variate Gaussian distribution, and 0 = {tt, ..., fip, Zi,..., Zl}. We cluster 
y^,..., into L clusters using K-Means clustering and randomly split the members of every cluster into k 
subsets. This ensures that full-data are split into k subsets such that the mixture proportions are represented 
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in every subset. Let Uji, • • •, yjra represent the data on subset j (j = 1,..., k). The hierarchical model for 
the data on subset j is 

L 

I Zji,0 - (44) 

h=l 

7t ~ Dirichlet(L“^,..., L“^), | ~ (0) IOOIh); Ik ~ Inverse-Wishart(2,4Ip), h = l,...,L, 

where 2 is the prior degrees of freedom and 4Ip is the scale matrix of the Inverse-Wishart distribution. 

The posterior distribution of 0 after stochastic approximation is derived using standard arguments for 
finite mixture of Gaussians. The likelihood given ..., and latent variables Zji,..., Zjra is 


h=l 


(45) 


where l(zji = h) is 1 if zji = h and 0 otherwise and (Jhj = l(zji = h.). For stochastic approxima¬ 

tion, we raise Lj in ( [43] ) to the power y and obtain 


h=l 


(46) 


If we use as the likelihood in ( |43] ), then the prior for 0 in ( |43l ) and simple extensions of standard arguments 
for hnite mixture of Gaussians imply that the analytic form of full conditional densities of the parameters 
are as follows. Dehne 


hj ={i: l(zji =H) = 1, i = l,...,m}, hhj = ^ ^ hji (H = 1,..., L), 

^ ’ ieHj 

and the complete conditionals of the parameters are 

/ m m \ 

7tj I rest ~ Dirichlet ( l(zji = 1) + ... ,y ^ l(zji = L) + j 

Pj I rsst ~ Normal 

IjH I rest ~ Inverse-Wishart I ytJhj -h 3, ^ (yjt -yn^ )(yii -hUj + 


(47) 




\ 0.01 TySlij^’^L 0.01+ytjhj 


i—1 

ih 


O.Ol-ytthj^ 


iehi 


0.01+ySHj 


ly+yhj +4 It 


Zji I rest ~ L Pjh^H) Pj 


iH = 


7tiH7^p(yji I PjHi^jh) 


(48) 


h=l ^H=l I l^jh’^jh,^ 

for h, = 1,..., L and i = 1,..., m; see Chapter 9 in [Bishop] ( |2006 l for details. 

Ah full conditionals are analytically tractable in terms of standard distributions. The Gibbs sampler 
iterates between the following four steps: 


1. Sample Ttj from Dirichlet distribution in ( [48] ). 

2. Sample Pjp^ from Normal distribution in ( [43] ) for h. = 1,..., L. 


3. Sample Zjh from Inverse-Wishart distribution in (|48|) for h = 1,..., L. 


4. Sample zji from categorical distribution in ( |48) for i = 1,..., m. 
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2.2 Simulated data analysis: Linear mixed effects model 

The conditional densities of |3 and L in two steps. First, the sampling model of the linear mixed effects 
model implies that the likelihood of ith observation in jth subset is 


l-ii = 


^ni(yji I Xi |3+ZiUt,T^InJ3Sfq(Ui | 0,I)dUi = 3srni(yji I Xi 


^,ZiLZl + xHr 


(49) 


This implies that likelihood of (3 and L after stochastic approximation is 

m m 

= Yl = n I p- . ( 50 ) 

1=1 1=1 

Second, the jth subset posterior distribution of |3 and L after stochastic approximation is calculated using 
LT as the likelihood and the same priors for |3 and L as in the sampling model of the linear mixed effects 


model. Instead of finding the analytic form of t 

le posterior density, we use the increment_log_prob 

function in Stan ( 

Stan Development Team 

2014 

1 to specify that the likelihood of as ( 

50 1 and use 


the default priors of |3 and L in Stan to obtain samples of (|3,1) from the jth subset posterior after stochastic 
approximation. Table ( fTO] ) describes the accuracies of CMC, SDP, VB, and WASP in approximating the full 
data marginal posterior distributions of fixed effects in the simulation. 


Table 10: Accuracies of the approximate posteriors for (3 in linear mixed effects model simulation. The 
accuracies are averaged over 10 simulation replications and across all elements of p. Monte Carlo errors are 
in parenthesis. VB, variational Bayes; CMC, consensus Monte Carlo; SDP, semiparametric density product; 
WASP, Wasserstein posterior. 



p=4 

p =80 

VB 

0.96 (0.01) 

0.96 (0.01) 


k= 10 

7^ 

II 

to 

o 

k= 10 

k = 20 

CMC 

0.96 (0.01) 

0.95 (0.02) 

0.95 (0.02) 

0.94 (0.03) 

SDP 

0.96 (0.01) 

0.95 (0.02) 

0.90 (0.06) 

0.90 (0.06) 

WASP 

0.95 (0.01) 

0.96 (0.02) 

0.95 (0.02) 

0.94 (0.03) 


2.3 Simulated data analysis: Probablistic parafac model 

The derivation of modified full conditional densifies of unknown parameters involves one key modification 


in fhe Gibbs sampling algorifhm of Dunson and Xing (20091. Assuming Zji is given, fhe confribufion of ifh 


observation in jfh subsef fo fhe likelihood after stochastic approximafion is 


L 




(D- 

H 


t* 

H=l) 


■, (iph ■ • •, =I n n n^ht 


l(x„ =l,z„=h) 


, h=l 


q=lt=l 


where l(xjiq = I, Zji = h) is 1 if bofh conditions are frue and is 0 ofherwise and I* is fhe maximum number 
of afoms in fhe sfick breaking represenfafion for fhe disfribufion of Zji. The condifional posferior densify of 
affer sfochasfic approximafion is proportional fo 

t=i 1=1t=i t=i 


1=1^) 




(q)aji+YLr;i^i(xjiq=i,Zji=h)-i 

ht > 
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which implies that 


I rest ~ Dirichlet (aqi +yl(xjiq = l,Zji = h),..., aqdq +yl(xjtq = dq,Zji = h,)) 


(51) 


for q = 1,..., p and h. = 1,..., I*. The conditional densities of remaining parameters follow from Section 
3.1 of Dunson and Xing ( |2009 1: 


VjH I rest ~ Beta(l + y ^ l(zji = h), a + y ^ l(zji > h,)), 

i=l i=l 

t* 

oij I rest ~ Gamma(aa + l*,ba — ^ log(l — Vjh))- 


(52) 

(53) 


h=l 


Finally, we update the posterior density of responsibility of every observation as 


Zji I rest■ 




h) Pjti — 




V] -V . . 

I LA-j-^q 


H=1 




(h, = 1,... ,r; i = 1,... ,m), (54) 


where WjH = Vjh Ot<H(l ~^jt)- The conditional posterior densities without stochastic approximation are 
obtained by substituting y = 1 and m = n in the full conditionals ( [51] ) - (531. 

All full conditionals are analytically tractable in terms of standard distributions. The Gibbs sampler 
iterates between the following four steps: 

1. Sample from Dirichlet distribution in ( [5T] ) for q = 1,... ,p and h = 1,..., I*. 


2. Sample Vjh from Beta distribution in ([52|) for h = 1,..., I*. 


3. Sample oij from Gamma distribution in ([53]). 


4. Sample zji from categorical distribution in ( |54| ) for i = 1,..., m. 
We fix Oq!, = 1/dq for q = 1,... ,p and I = 1,..., dq. 


2.4 Real data analysis: MovieLens data 


Table 11 describes the accuracies of CMC, SDP, VB, and WASP in approximating the full data marginal 
posterior distributions of fixed effecfs in the MovieLens data analysis. 


Table 11: Accuracies of the approximate posteriors of the fixed effecfs in the linear mixed effects model for 
MovieLens data. The accuracies are averaged over 10 replications. Monte Carlo errors are in parenthesis. 
CMC, consensus Monte Carlo; SDP, semiparametric density product; WASP, Wasserstein posterior. 



(^Action 

^Children — Action 

P Comedy — Action 

3 Drama — Action 

(^Popularity 

(^Previous 

CMC 

0.95 (0.02) 

0.93 (0.02) 

0.92 (0.02) 

0.95 (0.02) 

0.94 (0.02) 

0.95 (0.03) 

SDP 

0.93 (0.03) 

0.92 (0.03) 

0.92 (0.04) 

0.93 (0.04) 

0.92 (0.04) 

0.94 (0.02) 

VB 

0.00 (0.00) 

0.34 (0.05) 

0.64 (0.05) 

0.15(0.02) 

0.01 (0.00) 

0.00 (0.00) 

WASP 

0.96 (0.01) 

0.95 (0.01) 

0.95 (0.02) 

0.96 (0.01) 

0.96 (0.01) 

0.96 (0.01) 
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3 MATLAB Code for solving the WASP linear program 


For ease of illustration, we eonsider the problem of estimating an atomie approximation of the WASP 
using three subset posterior distributions eonsisting of 50, 75, and 100 posterior samples from Gaus¬ 
sian distributions with means = (3,2)^, ^2 = (2,3)^, and [ 0.3 = (3,3)^ and eovarianee matriees 
Zi = Z 2 = 2^3 = Z, where On = 1, ^22 = 3, and a 2 i = 1.5. The MATLAB funetion WASP, whieh esti¬ 
mates the atomie approximation of the WASP, requires two inputs. First, the support of the WASP. Seeond, 
the distanee matrix between the atoms of the WASP and every subset posterior distribution. We assume that 
the WASP is supported on a grid of atoms estimated from all subset posterior samples, where mesh-size 
of the grid is speeified by grdsize. We use the grid of atoms of the WASP to obtain the three distanee 
matriees between the atoms of the WASP and every subset posterior distribution. The MATLAB funetion 
WASP uses the grid of atoms and distanee matriees as inputs and estimates the weights of every atom in the 
support of the WASP. The empirieal measure obtained using the atoms and their estimated weights is an 
atomie approximation of the true analytieally intraetable WASP. 

rand('seed',0); 

% Size of the grid 
grdsize = 50; 

% var-covar matrix for Multivariate Normal Distribution 
sig = [1 1.5; 1.5 3]; 

% means for 3 subset posteriors 
mul = [3 2]; 
mu2 = [23]; 
mu3 = [3 3]; 

% atoms for 3 subset posteriors; they only differ in their means, 
spostfl} = mvnrnd(mul, sig, 50); 

spost{2} = mvnrnd(mu2, sig, 75); 

spost{3} = mvnrnd(mu3, sig, 100); 

% atoms for the WASP by forming a grid 

Ibdl = min{cellfun(@(x) x(l), cellfun(@(x) min{x), spost,'UniformOutput', false))); 

lbd2 = min{cellfun(@(x) x(2), cellfun(@(x) min(x), spost,'UniformOutput', false))); 

ubdl = max{cellfun(@(x) x(l), cellfun(@(x) max{x), spost,'UniformOutput', false))); 

ubd2 = max{cellfun(@(x) x(2), cellfun(@(x) max{x), spost,'UniformOutput', false))); 

[opostx, oposty] = meshgrid{linspace(Ibdl, ubdl, grdsize), linspace{lbd2, ubd2, grdsize)); 
opost = [opostx{:) oposty (:)]; 

% calculate the pair-wise sq. euclidean distance between the atoms of subset 

% posteriors and BarPost atoms 

mil = diag{spost{1} * spost{l}'); 

m22 = diag{spost{2} * spost{2}'); 

m33 = diag{spost{3} * spost{3}'); 

mOO = diag{opost * opost'); 

mOl = opost * spost{1}'; 
m02 = opost * spost{2}'; 
m03 = opost * spost{3}'; 

% calculate distance between atoms 

dOl = bsxfun{@plus, bsxfun(@plus, -2 * mOl, mil'), mOO); 

d02 = bsxfun{@plus, bsxfun (@plus, -2 * m02, m22'), mOO); 

d03 = bsxfun{@plus, bsxfun(@plus, -2 * m03, m33'), mOO); 

% initialize the wts b_l ... b_K for subset posteriors atoms; see Eq. 21 

b = cellfun{@{x) ones(size(x, 1), 1) / size(x, 1), spost, 'UniformOutput', false); 

colMeasure = b; 
distMat{l} = dOl; 
distMat{2} = d02; 
distMat{3} = d03; 

tic 

[optSol, optObj, exitflag, output, lambda] = WASP{colMeasure, distMat); 
toe 

function [optSol, obj, exitflag, output, lambda] = WASP(colMeasure, distMat) 

%% WASP function to calculate the Wasserstein Barycenter of a list of 
% empirical measures and pairwise distance between atoms. 

% Input: 

% ====== 

% colMeasure: 

% cell containing the wts of atoms of empirical measures. 

% 

% distMat: 

% cell containing pairwise distance between atoms of empirical measures and 
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% the Wasserstein barycenter. 

% 

% Output: 

% ====== 

% output of solving the LP (31) in the manuscript using 'linprog' 
nsubset = length(distMat); 

vecColMeasure = cell2mat(cellfun(@(x) x', colMeasure, 'UniformOutput', false))'; 
nsample = size(vecColMeasure, 1); 

vecDistMat = cell2mat(distMat); 
vecDistMat = vecDistMat(:); 

fmatCell = cellfun(@(x) kron(ones(l, size(x, 1)), eye(nsample)), colMeasure, 'UniformOutput', false); 
fmat = blkdiag(fmatCell{:}); % F 

hmatCell = cellfun(@(x) kron(eye(size(x, 1)), ones(l, nsample)), colMeasure, 'UniformOutput', false); 
hmat = blkdiag(hmatCell{:}); % H 

gmat = kron(ones(nsubset, 1), eye(nsample)); % G 

% Aeq in the matlab linprog function 
aMat = [zeros (1, nsample''2) ones(l, nsample); 
fmat -gmat; 

hmat zeros(nsample, nsample) 

]; 

% beq in the matlab lingprog function 
bVec = [1; 

zeros(nsubset * nsample, 1); 
vecColMeasure 
]; 

% f in the matlab linprog function 
costVec = [vecDistMat; 

zeros(nsample, 1)]; 

% upper bds = 1 and lower bds = 0 
Ibd = zeros (nsample''2 + nsample, 1); 
ubd = ones (nsample''2 + nsample, 1); 

[optSol, obj, exitflag, output, lambda] = linprog(costVec, [], [], aMat, bVec, Ibd, ubd); 
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